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ABSTRACT 

We use multi-scale smoothed particle hydrodynamic simulations to study the inflow of gas 
from galactic scales (<~ lOkpc) down to < 0.1 pc, at which point the gas begins to resemble 
a traditional, Keplerian accretion disk. The key ingredients of the simulations are gas, stars, 
black holes (BHs), self-gravity, star formation, and stellar feedback (via a subgrid model); 
BH feedback is not included. We use <~ 100 simulations to survey a large parameter space of 
galaxy properties and subgrid models for the interstellar medium physics. We generate initial 
conditions for our simulations of galactic nuclei (< 300 pc) using galaxy scale simulations, 
including both major galaxy mergers and isolated bar-(un)stable disk galaxies. For sufficiently 
gas-rich, disk-dominated systems, we find that a series of gravitational instabilities generates 
large accretion rates of up to <~ 1 — lOMQyr -1 onto the BH (i.e., at < 0.1 pc); this is com- 
parable to what is needed to fuel the most luminous quasars. The BH accretion rate is highly 
time variable for a given set of conditions in the galaxy at <~ kpc. At radii > 10 pc, our sim- 
ulations resemble the "bars within bars" model of Shlosman et al, but we show that the gas 
can have a diverse array of morphologies, including spirals, rings, clumps, and bars; the duty 
cycle of these features is modest, complicating attempts to correlate BH accretion with the 
morphology of gas in galactic nuclei. At ~ 1 — 10 pc, the gravitational potential becomes 
dominated by the BH and bar-like modes are no longer present. However, we show that the 
gas can become unstable to a standing, eccentric disk or a single-armed spiral mode (m = 1), 
in which the stars and gas precess at different rates, driving the gas to sub-pc scales (again 
for sufficiently gas-rich, disk-dominated systems). A proper treatment of this mode requires 
including star formation and the self-gravity of both the stars and gas (which has not been 
the case in many previous calculations). Our simulations predict a correlation between the 
BH accretion rate and the star formation rate at different galactic radii. We find that nuclear 
star formation is more tightly coupled to AGN activity than the global star formation rate of a 
galaxy, but a reasonable correlation remains even for the latter. 
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1 INTRODUCTION 

The inflow of gas into the central parts of galaxies plays a crit- 
ical role in galaxy formation, ultimately generating phenomena 
as diverse as bulges and spheroidal galaxies, starbursts and ultra- 
luminous infrared galaxies (ULIRGs), nuclear stellar clusters, and 
accretion onto super-massive black holes (BHs). The discovery, in 
the past decade, of tight correlations between black hole mass and 
host spheroid properties including mass (Kormendy & Richstone 
1995; Magorrian et al. 1998), velocity dispersion (Ferrarese & Mer- 
ritt 2000; Gebhardt et al. 2000), and binding energy or potential 
well depth (Hopkins et al. 2007; Aller & Richstone 2007) implies 
that these phenomena are tightly coupled. 

It has long been realized that bright, high-Eddington ratio ac- 
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cretion (i.e., a quasar) dominates the accumulation of mass in the 
supermassive BH population (Soltan 1982; Salucci et al. 1999; 
Shankar et al. 2004). In order to explain the existence of black holes 
with masses ~ 1O 9 M , the amount of gas required is comparable 
to that contained in entire large galaxies. Given the short lifetime of 
the quasar phase < 10 s yr (Martini 2004), the processes of interest 
must deliver a galaxy's worth of gas to the inner regions of a galaxy 
on a relatively short timescale. 

There is also compelling evidence that quasar activity is pre- 
ceded and/or accompanied by a period of intense star formation in 
galactic nuclei (Sanders et al. 1988a,b; Dasyra et al. 2007; Kauff- 
mann et al. 2003). The observed properties of bulges at z ~ inde- 
pendently require that dissipative processes (gas inflow) dominate 
the formation and structure of the inner ~kpc (Ostriker 1980; Carl- 
berg 1986; Gunn 1987; Kormendy 1989; Hernquist et al. 1993). 
Hopkins et al. (2009a,d, 2008a) showed that this inner dissipational 
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component can constitute a large fraction ~ 5 — 30% of the galaxy's 
mass, with stellar (and at some point probably gas) surface densi- 
ties reaching ~ 10 n ~' 2 MQkpc~ 2 . 

On large (galactic) scales, several viable processes for initi- 
ating such inflows are well-known. Major galaxy-galaxy mergers 
produce strong non-axisymmetric disturbances of the constituent 
galaxies; such disturbances may also be produced in some mi- 
nor mergers and/or globally self-gravitating isolated galactic disks. 
Observationally, major mergers are associated with enhancements 
in star formation in ULIRGs, sub-millimeter galaxies, and pairs 
more generally (e.g. Sanders & Mirabel 1996; Schweizer 1998; Jo- 
gee 2006; Dasyra et al. 2006; Woods et al. 2006; Veilleux et al. 
2009). Numerical simulations of mergers have shown that when 
such events occur in gas-rich galaxies, resonant tidal torques lead to 
rapid inflow of gas into the central ~kpc (Hernquist 1989; Barnes 
& Hernquist 1991, 1996). The resulting high gas densities trigger 
starbursts (Mihos & Hernquist 1994, 1996), and are presumed to 
feed rapid black hole growth. Feedback from the starburst and a 
central active galactic nucleus (AGN) may also be important, both 
for regulating the BH's growth (Di Matteo et al. 2005; Hopkins 
et al. 2005; DeBuhr et al. 2009; Johansson et al. 2009b) and for 
shutting down future star formation (Springel et al. 2005a; Johans- 
son et al. 2009a; see, however, DeBuhr et al. 2009). 

However, the physics of how gas is transported from ~ 1 kpc 
to much smaller scales remains uncertain (e.g., Goodman 2003). 
Typically, once gas reaches sub-kpc scales, the large-scale torques 
produced by a merger and/or large-scale bar/spiral become less ef- 
ficient. In the case of stellar bars or spiral waves, there can even be 
a "hard" barrier to further inflow in the form of an inner Linblad 
resonance, if the system has a non-trivial bulge. In mergers, the co- 
alescence of the two systems generates perturbations on all scales, 
and so allows gas to move through the resonances, but the pertur- 
bations relax rapidly on small scales, often before gas can inflow. 

Local viscous stresses - which are believed to dominate angu- 
lar momentum transport near the central BH (e.g., Balbus & Haw- 
ley 1998) - are inefficient at radii > 0.01 — 0.1 pc (e.g., Shlosman 
& Begelman 1989; Goodman 2003; Thompson et al. 2005). It is 
in principle possible that some molecular clouds could be scattered 
onto very low angular momentum orbits, but even the optimistic 
fueling rates from this process are generally insufficient to produce 
luminous quasars (see e.g. Hopkins & Hernquist 2006; Kawakatu 
& Wada 2008; Nayakshin & King 2007). As a consequence, many 
models invoke some form of gravitational torques ("bars within 
bars"; Shlosman et al. 1989) to continue transport to smaller radii. 
As gas is driven into the central kpc by large-scale torques, it will 
cool rapidly into a disky structure; if this gas reservoir is massive 
enough, the gas will be self-gravitating and thus again vulnerable 
to global instabilities (e.g., the well-known bar and/or spiral wave 
instabilities) that can drive some of the gas to yet smaller radii. 

To date, numerical simulations have seen the formation of 
such secondary bars in some circumstances, such as in adaptive 
mesh refinement (AMR) simulations of galaxy formation (Wise 
et al. 2007; Levine et al. 2008; Escala 2007) or particle-splitting 
smoothed particle hydrodynamics (SPH) simulations of some ide- 
alized systems (Escala et al. 2004; Mayer et al. 2007). These stud- 
ies have served as a critical "proof of concept." However, these ex- 
amples have generally been limited by computational expense to 
studying a single system at one instant in its evolution, and thus it 
is difficult to assess how the sub-pc dynamics depends on the large 
parameter space of possible inflow conditions from large radii and 
galaxy structural parameters. 

Alternatively, some simulations simply take an assumed 



small-scale structure and/or fixed inflow rate as an initial/boundary 
condition, and study the resulting gas dynamics at small radii (e.g. 
Schartmann et al. 2009; Dotti et al. 2009; Wada & Norman 2002; 
Wada et al. 2009). These studies have greatly informed our under- 
standing of nuclear obscuration on small scales (the "torus"), and 
the role of stellar feedback in determining the structure and dy- 
namics of the gas at these radii; it is, however, unclear how to relate 
this small-scale dynamics to the larger-scale properties of the host 
galaxy. This is critical for understanding black hole growth and nu- 
clear star formation in the broader context of galaxy formation. 

Observationally, a long standing puzzle has been that many 
systems, especially those with weaker inflows on large scales (e.g. 
bar or spiral wave-unstable disks with some bulge, as opposed to 
major mergers), exhibit no secondary instabilities at ~ 0.1 — 1 kpc 
- in several cases, torques clearly reverse sign inside these radii 
(Block et al. 2001; Garcfa-Burillo et al. 2005). Whether this is 
generic, or the consequence of a low duty cycle, or the result of 
the large-scale inflows simply being too weak in these cases, is not 
clear. Moreover, even among systems that do show nuclear asym- 
metries, and that clearly exhibit enhanced star formation and lu- 
minous AGN, the observed features at smaller radii are very often 
not traditional bars. Rather, they exhibit a diverse morphology, with 
spirals quite common, along with nuclear rings, barred rings, occa- 
sional one or three-armed modes, and some clumpy/irregular struc- 
tures (Martini & Pogge 1999; Peletier et al. 1999; Knapen et al. 
2000; Laine et al. 2002; Knapen et al. 2002; Greene et al. 2009). 

Even if secondary bars or spirals are present at intermediate 
radii ~ 10 — 100 pc, it has long been recognized that they will 
cease to be important at yet smaller scales, when the potential be- 
comes quasi-Keplerian and the global self-gravity of the gas less 
important; this occurs as one approaches the BH radius of influ- 
ence, which is ~ lOpc in typical ~ L* galaxies (Athanassoula et al. 
1983, 2005; Shlosman et al. 1989; Heller et al. 2001; Begelman & 
Shlosman 2009). Indeed, in previous simulations and most analytic 
calculations, the "bars-within-bars" model appears to break down 
at these scales (see e.g. Jogee 2006, and references therein). How- 
ever, local angular momentum transport is still very inefficient at 
~ 10 pc, and the gas is still locally self-gravitating, and so should be 
able to form stars rapidly (e.g., Thompson et al. 2005). Understand- 
ing the physics of inflow through these last few pc, especially in a 
consistent model that connects to gas on galactic scales (~ lOkpc), 
remains one of the key open questions in our understanding of mas- 
sive BH growth. 

In this paper, we present a suite of multi-scale hydrodynamic 
simulations that follow gravitational torques and gas inflow from 
the kpc scales of galaxy-wide events through to < 0. 1 pc where the 
material begins to form a standard thin accretion disk. These sim- 
ulations include gas cooling, star formation, and self-gravity; feed- 
back from supernovae and stellar winds is crudely accounted for via 
a subgrid model. In order to isolate the physics of angular momen- 
tum transport, we do not include BH feedback in our calculations. 
We systematically survey a large range of galaxy properties (e.g., 
gas fraction and bulge to disk ratio) and gas thermodynamics, in 
order to understand how these influence the dynamics, inflow rates, 
and observational properties of gas on small scales in galactic nu- 
clei (~ 0.1 — 100 pc). Our focus in this paper is on the results of 
most observational interest: what absolute inflow rates, star forma- 
tion rates, and gas/stellar surface density profiles result from sec- 
ondary gravitational instabilities? What is their effective duty cy- 
cle? And what range of observational morphologies are predicted? 
In a future paper (Paper II) we will present a more detailed compar- 
ison between our numerical results and analytic models of inflow 
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and angular momentum transport induced by non-axisymmetric in- 
stabilities in galactic nuclei. 

The remainder of this paper is organized as follows. In § 2 we 
describe our simulation methodology, which consists of two levels 
of "re-simulations" using initial conditions motivated by galactic- 
scale simulations. In § 3 we present an overview of our results and 
show how a series of gravitational processes leads to gas transport 
from galactic scales to sub-pc scales. In §4 we quantify the resulting 
inflow rates and gas properties as a function of time and radius in 
the simulations. §5 summarizes the conditions required for global 
gravitational instability and significant gas inflow. In § 6 we show 
how the physics of accretion induced by gravitational instabilities 
leads to a correlation (with significant scatter) between star for- 
mation at different radii and BH accretion; we also compare these 
results to observations. In § 7 we summarize our results and discuss 
a number of their implications and several additional observational 
tests. Further numerical details and tests of our methodology are 
discussed in § A. In § B we show how the subgrid model of the 
ISM we use influences our results. 

2 METHODOLOGY 

We use a suite of hydrodynamic simulations to study the physics 
of gas inflow from ~ lOkpc to ~ 0.1 pc in galactic nuclei. In or- 
der to probe the very large range in spatial and mass scales, we 
carry out a series of "re-simulations." First, we simulate the dynam- 
ics on galaxy scales. Specifically, we use representative examples 
of gas-rich galaxy-galaxy merger simulations and isolated, moder- 
ately bar-unstable disk simulations. These are well-resolved down 
to ~ 100 — 500pc. We use the conditions at these radii (at several 
times) as the initial conditions for intermediate-scale re-simulations 
of the sub-kpc dynamics. In these re-simulations, the smaller vol- 
ume is simulated at higher resolution, allowing us to resolve the 
subsequent dynamics down to ~ lOpc scales - these re-simulations 
approximate the nearly instantaneous behavior of the gas on sub- 
kpc scales in response to the conditions at ~kpc set by galaxy-scale 
dynamics. We then repeat our re-simulation method to follow the 
dynamics down to sub-pc scales where the gas begins to form a 
standard accretion disk. 

Our re-simulations are not intended to provide an exact re- 
alization of the small-scale dynamics of the larger-scale simula- 
tion that motivated the initial conditions of each re-simulation (in 
the manner of particle-splitting or adaptive-mesh refinement tech- 
niques). Rather, our goal is to identify the dominant mechanism(s) 
of angular momentum transport in galactic nuclei and what param- 
eters they depend on. This approach clearly has limitations, espe- 
cially at the outer boundaries of the simulations; however, it also 
has a major advantage. By not requiring the conditions at small 
radii to be uniquely set by a larger-scale "parent" simulation, we 
can run a series of simulations with otherwise identical conditions 
(on that scale) but systematically vary one parameter (e.g., gas frac- 
tion or ISM model) over a large dynamic range. This allows us to 
identify the physics and galaxy properties that have the biggest ef- 
fect on gas inflow in galactic nuclei. As we will show, the diversity 
of behaviors seen in the simulations, and desire to marginalize over 
the uncertain ISM physics, makes such a parameter survey critical. 

This methodology is discussed in more detail below. First, we 
describe the physics in our simulations, in particular our treatment 
of gas cooling, star formation, and feedback from supernovae and 
young stars (§ 2.1). We then summarize the galaxy-scale simula- 
tions that are used to motivate the initial conditions for subsequent 
re-simulations (§ 2.2). The intermediate-scale re-simulations, and 



the methodology used to construct their initial conditions, are dis- 
cussed in § 2.3. Finally, we discuss the nuclear-scale resimulations, 
which are themselves motivated by the intermediate- scale resimu- 
lations (§ 2.4). 

2.1 Gas Physics, Star Formation, and Stellar Feedback 

The simulations were performed with the parallel TreeSPH code 
GADGET- 3 (Springel 2005), based on a fully conservative formula- 
tion of smoothed particle hydrodynamics (SPH), which conserves 
energy and entropy simultaneously even when smoothing lengths 
evolve adaptively (see e.g., Springel & Hernquist 2002; Hernquist 
1993; O'Shea et al. 2005). The detailed numerical methodology 
is described in Springel (2005), Springel & Hernquist (2003), and 
Springel et al. (2005b). 

The simulations include supermassive black holes (BHs) as 
additional collisionless particles at the centers of all progenitor 
galaxies. In our calculations the BH's only dynamical role is via 
its gravitational influence on the smallest scales ~ 1 — 10 pc. In 
particular, to cleanly isolate the physics of gas inflow, we do not in- 
clude the subgrid models for BH accretion and feedback that have 
been used in previous works (e.g., Springel et al. 2005b). During 
a galaxy merger, the BHs in each galactic nucleus are assumed to 
coalesce and form a single BH at the center of mass of the system 
once they are within a single SPH smoothing length of one another 
and are moving at a relative speed lower than both the local gas 
sound speed and relative escape velocities. 

In our models, stars form from the gas using a prescription 
motivated by the observed Kennicutt (1998) relation. Specifically, 
we use a star formation rate per unit volume p\ oc p 3//2 with the 
normalization chosen so that a Milky-way like galaxy has a total 
star formation rate of about 1 Mq yr~ 1 . 

The precise slope, normalization, and scatter of the Schmidt- 
Kennicutt relation, and even whether or not such a relation is gen- 
erally applicable, are somewhat uncertain on the smallest spatial 
scales we model here. This is especially true when the dynamical 
times become short relative to the main-sequence stellar lifetime 
(fdyn ~ 10 5 — 10 6 yr in the smallest regions simulated). Nonethe- 
less, there is some observational and physical motivation for the 
"standard" parameters we have adopted, even at high surface den- 
sities. For the densest star forming galaxies, observational studies 
favor a logarithmic slope ~ 1 .7 for the relation between E» and E g 
(Bouche et al. 2007), not that different from what our model imple- 
ments. In addition, Tan et al. (2006) and Krumholz & Tan (2007) 
show that local observations imply a constant star formation effi- 
ciency in units of the dynamical time (i.e. p„ oc p 1 ' 5 ) at all densities 
observed, n ~ 10'~ 6 cm~ 3 - the highest gas densities in these stud- 
ies are comparable to the highest gas densities in our simulations 
(~ 10 s Mq of gas inside ~ 10 pc). Finally, Davies et al. (2007) & 
Hicks et al. (2009) estimate the star formation rate (SFR) and gas 
surface densities in AGN on exactly the small scales of interest here 
(~ 1 — lOpc); they find a SFR-density relation continuous with that 
implied at "normal" galaxy densities. 

To understand the possible impact of uncertainties in the 
Schmidt-Kennicutt relation on our conclusions, we have adjusted 
the slope d In p\ /d In p adopted in our simulations between 1.0 — 2.0 
in a small set of test runs, fixing the star formation rate at MW- 
like surface densities of ~ 10 9 MQkpc~ 2 where the observational 
constraints are tight. This amounts to varying the absolute star for- 
mation efficiency on the smallest resolved scales by a factor of 
> 100; qualitatively, this could presumably mimic a wide variety 
of different physics associated with stellar feedback and star for- 
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Figure 1. Effective equation of state of the ISM in our simulations; this accounts for the effects of stellar feedback that are not resolved in our calculations. We 
plot the effective sound speed c s (i.e., the turbulent speed generated by feedback) versus average ISM density n. The q cos = case is an isothermal "floor" at 
c s = 10km s _ 1 . The q eos = 1 line is the "maximal feedback" model in Springel & Hernquist (2003), in which the ISM is multi-phase above a minimum n and 
all supernova energy goes into pressurizing the "diffuse" medium. Intermediate q ms interpolate between the two (eq. 1). For each q sas , we show an analytic 
curve for the equilibrium c s (n), and simulation results (colored points). We compile measurements of the turbulent or non-thermal velocities in observed 
systems (black points): local ULIRGs (Downes & Solomon 1998, circles; large/small for the inner/outer observed radii in each), the nuclei of merging LIRGs 
(Bryant & Scoville 1999, star), the core of M82 (Westmoquette et al. 2007, x), the central 400 pc of NGC 6240 (Tacconi et a!. 1999, triangle), nuclear 
clumps in NGC 6240 (Iono et al. 2007, pentagon), and the maser disk of NGC 3079 (Kondratko et al. 2005, small squares). The canonical spiral disk velocity 
dispersion is also shown (large square). At high-redshifts, we show recent observations of "normal" star-forming disks from Forster Schreiber et al. (2006, 
z ~ 2; diamond) and Lemoine-Busserolle et al. (2009, z = 3 — 4; asterisk), and luminous sub-millimeter galaxies from Tacconi et al. (2006, inverted triangle). 
The observations favor a median q tos ~ 0. 125 — 0.3, which we adopt as our "fiducial" models. 



mation. This variation can, unsurprisingly, have a dramatic affect 
on the quasi-equilibrium gas densities at small radii, which are set 
by gas inflow balancing star formation. However, even over this 
large range of star formation efficiencies, the qualitative behavior 
of the angular momentum transport and gas inflow does not change 
dramatically; the gas dynamics in a low-star formation efficiency 
run is similar to that in a run with much higher initial gas content 
but also higher star formation efficiency. As a result, although the 
absolute star formation efficiency is clearly somewhat uncertain, 
we do not believe that this qualitatively affects our conclusions. 

The largest uncertainties in our modeling stem from the treat- 
ment of the interstellar medium (ISM) gas physics and the im- 
pact of stellar feedback on the ISM. Our simulations are rela- 
tively coarse and average over many star-forming clumps, HII re- 
gions, supernova remnants, etc. As a result, the simulations use a 
sub-resolution model of a multiphase ISM in which the gas has a 
sound speed much larger than its true thermal velocity (Springel 
& Hernquist 2003). Our assumption is that the large-scale grav- 
itational torques produced by bars, spiral waves, and other non- 



axisymmetric features, will not depend critically on the small-scale 
structure of the ISM; although we believe that this is qualitatively 
correct, more detailed calculations will be required to ultimately 
assess this assumption. The key role of stellar feedback in this 
model is to suppress the runaway fragmentation and clumping of 
gas on small scales. In reality, this likely occurs via turbulence gen- 
erated by stellar feedback and via the disruption of star clusters and 
molecular clouds (e.g., Murray et al. 2009). In our model, all of 
this physics is "accounted for" by the large effective sound speed, 
which increases the Jeans and Toomre masses, thus suppressing the 
formation of small-scale structure. 

Figure 1 shows the range of effective sound speeds c s in our 
calculations as a function of the ISM density n, compared to a num- 
ber of observational constraints (large symbols). The solid lines in 
Figure 1 represent analytic approximations to the equation of state, 
while the small colored points are representative results from sim- 
ulations that also include adiabatic cooling/heating and shock heat- 
ing. We adopt a parameterization of the sound speed as a function 
of density - i.e. the effective equation of state for the ISM gas - 



© 0000 RAS, MNRAS 000, 000-000 



How Do Massive Black Holes Get Their Gas ? 5 



following Springel et al. (2005b); Springel & Hernquist (2005); 
Robertson et al. (2006a,b). With this model, we can interpolate 
freely between two extremes using a parameter q eofi . At one ex- 
treme, the gas has an effective sound speed of 10kms~', motivated 
by, e.g., the observed turbulent velocity in atomic gas in nearby 
spirals or the sound speed of low density photo-ionized gas; this 
is the "no-feedback" case with q eos = 0. The opposite extreme, 
<7eos = 1, represents the "maximal feedback" sub-resolution model 
of Springel et al. (2005b), based on the multiphase ISM model 
of McKee & Ostriker (1977); in this case, 100% of the energy 
from supernovae is assumed to stir up the ISM. This equation of 
state is substantially stiffer, with effective sound speeds as high as 
~ 200km s _1 . Note that at the highest densities, c„ begins to de- 
cline in all of the models (albeit slowly), as the efficiency of star 
formation asymptotes but cooling rates continue to increase. 

Varying q eos between these two extremes amounts to varying 
the effective sound speed of the ISM, with the interpolation 

Cs = \J <7eos c? [q = 1 ] + ( 1 - <7eo S ) c? [q = 0] . (1) 

The resulting sound speeds for q eos — 0.125 and 0.25 are shown in 
Figure 1; these correspond to more moderate values of c s ~ 30 — 
100km s _1 for the densities of interest. 

Figure 1 compares these models to observations of the turbu- 
lent (non-thermal) velocities in atomic and molecular gas in a num- 
ber of systems (large symbols). At low mean densities, n < 0.3 — 
lcm -3 , the turbulent velocity in nearby spirals is ~ lOkrns -1 . 
Downes & Solomon (1998) present a detailed study of a number 
of luminous local starbursts that have significantly higher mean 
densities; they decompose the observed molecular line profiles into 
bulk (e.g., rotation) and turbulent motions. We plot their determi- 
nation of the mean density and turbulent velocities in each system 
at several radii. We also show the results of similar observations 
of the core of M82 (Westmoquette et al. 2007), additional nearby 
luminous infrared galaxies (Bryant & Scoville 1999), NGC 6240 
(Tacconi et al. 1999; Iono et al. 2007), and luminous starbursts at 
high redshift, z ~ 2 — 3 (Forster Schreiber et al. 2006; Tacconi et al. 
2006; Lemoine-Busserolle et al. 2009); at the highest densities, we 
also show the random velocities observed in the nuclear maser disk 
in the nearby Seyfert 2 galaxy NGC3079 (Kondratko et al. 2005). 

The observational results in Figure 1 favor models with q eos ~ 
0.1 —0.3, albeit with significant scatter. We thus take these values 
of q e0 s as our "standard" choices, although we have carried out nu- 
merical experiments over the entire range q eos = — 1 . Note that the 
observations clearly do not support a simple no-feedback (g eos = 0) 
model. Within the range q sos ~ 0. 1 — 0.3, our results on AGN fuel- 
ing are not particularly sensitive to the precise value of q t0!} . More- 
over, the functional form c s (n) is also not crucial: simulations using 
a constant c s ~ 50kms~' yield similar results. However, our simu- 
lations with q eos = and g eos = 1 predict results that are inconsistent 
with observations of galactic nuclei - thus, our results themselves 
favor q^s ^0.1 — 0.3 (see Appendix B). 

For the gas densities of interest in this paper, the precise form 
of the cooling law does not significantly affect our conclusions. 
This is because the cooling time is almost always much shorter 
than the local dynamical time (typical t coo i ~ 10~ 6 — 10~ 4 fd yn ). As 
a result the "sound speed" of the gas is nearly always pinned to 
the subgrid 'turbulent' value discussed above (this is why the nu- 
merical points in Fig. 1 are so close to the analytic models). This 
is true even when the gas is optically thick to the infrared radia- 
tion produced by dust, as can readily occur in the central ~ 100 
pc: the cooling time (diffusion time) is still much less than the dy- 
namical time in the optically thick limit for the radii that we re- 



solve (e.g., Thompson et al. 2005). We have, in fact, experimented 
with alternative cooling rate prescriptions: including or excluding 
metal-line cooling, uniformly increasing or decreasing the cooling 
rate by a factor of ~ 3, and in the most extreme case, assuming 
instantaneous gas cooling (any gas parcel above the cooling floor 
is assumed to immediately radiate the excess energy in a single 
timestep). We do not see any significant changes in our results with 
these variations, simply because the gas always cools rapidly in 
our calculations (in contrast, in regimes such as the a-disk where 
the cooling time is comparable to the dynamical time, the details 
of the cooling function can have a significant effect; see Gammie 
2001; Nayakshin et al. 2007; Cossins et al. 2009). If the effective 
minimum c s comes from turbulent velocities, then the OeffectiveO 
cooling time around this SSoor should be given by the turbulent de- 
cay time, which can be comparable to the dynamical time (Begel- 
man & Shlosman 2009); this is not included in our calculations. In 
the presence of such an effective cooling time, local gravitational 
instability may lead to tightly wound spirals as opposed to frag- 
mentation into star-forming clumps. These could be important for 
angular momentum transport at some radii. 

To conclude our discussion of the ISM physics in our sim- 
ulations, it is important to reiterate that the key role of the sub- 
resolution sound speed c s is that it determines the local Jeans and 
Toomre criteria, and thus the physical scale on which gravitational 
physics dominates. The Jeans mass for a disk of surface density 
£ and sound speed c s is given by Mj = (tyc 4 s )/(4G 2 £). For the 
outer regions of a galactic disk with S ~ 10 s — 10 9 Mq kpc~ 2 and 
c s ~ 10km s _1 , Mj ~ 10 6 Mq, comparable to that of a molecular 
cloud; the corresponding Jeans length is tens of pc, comparable to 
that of massive molecular cloud complexes in galaxies. Thus our 
sub-resolution model is effectively averaging over discrete molec- 
ular clouds and star clusters in galaxies. Large-scale inflows can 
increase the surface density in the central regions of galaxies, but 
c s also rises. In our models with q eos ~ 0. 1 — 0.3, the Jeans mass re- 
mains roughly similar down to ~pc scales, but as a result the Jeans 
length is significantly smaller in galactic nuclei where the ambient 
density is much higher. These physical mass and size-scales mo- 
tivate the numerical resolution in our simulations; in all cases, we 
ensure that the resolution is sufficient to formally resolve the Jeans 
mass and length. Higher resolution simulations may be numerically 
achievable, but can provide only minimal gains in the "reality" of 
the simulation without a corresponding increase in the sophistica- 
tion of the ISM model. 

2.2 Large Scale Galaxy Mergers and Bars: 100 kpc to 100 pc 

Our galaxy-scale simulations motivate the initial conditions chosen 
for the smaller-scale re-simulation calculations described in §2.3 & 
2.4. The galaxy-scale simulations include isolated disks (both glob- 
ally stable and bar unstable) and galaxy-galaxy mergers. We will 
ultimately focus on a few representative examples, but we chose 
those having surveyed a large parameter space. These simulations 
and the methodology used for building the initial galaxies are de- 
scribed in more detail in a series of papers (see e.g. Di Matteo et al. 
2005; Robertson et al. 2006b; Cox et al. 2006; Younger et al. 2008). 
We briefly review the key points here. 

For each simulation, we generate one or two stable, isolated 
disk galaxies, each with an extended dark matter halo with a Hern- 
quist (1990) profile, an exponential disk of gas and stars, and an 
optional stellar bulge. The initial systems are chosen to be consis- 
tent with the observed baryonic Tully-Fisher relation and estimated 
halo-galaxy mass scaling laws (Bell & de Jong 2001; Kormendy & 
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Table 1. Galaxy-Scale Simulations 



Simulation 


e 




/gas h/R B/T(<R) S d (0) 


S b (0) 


M tot 


(< R) [M ] 


Name 


[pc] 




(500 pc) (500 pc) (300 pc) (lkpc) [M kpc" 


- 2 ] 


(300 pc) 


(lkpc) 


(1) 


(2) 


(3) 


(4) (5) (6) (7) 


(8) 




(9) 


Merger & Isolated Galaxy Simulations: Parameter Studies 




20,50,100,150 


0.0-1.0 


0.05-1.0 0.05-0.35 0.1-0.9 0.01-0.5 1.0e9-1.0ell 1.0e9-1.0el2 


1.0e8-1.0e 


11 1.0e8-1.0ell 


Typical Gas-Rich ~ L* Merger: Initial Conditions 


b3ex(ic) 


10", 50 


0.25 


0.80 0.25 0.8 0.4 2.5e8 


2.5e8 


1.4e8 


1.9e9 


Typical Gas-Rich ~ L* Merger: At Coalescence 


b3ex(co) 


10", 50 


0.25 


0.45 0.33 0.15 0.25 2.0el0 


2.6e9 


4.4e9 


2.5el0 


Typical Gas-Rich ~ L» Merger: 10 s yr Post-Coalescence 


b3ex(po) 


10", 50 


0.25 


0.08 0.37 0.10 0.15 2.3ell 


9.1e9 


2.8el() 


4.6el0 


Bar-Unstable ~ L„ Disk: Initial Conditions 


barex(ic) 


10 a ,50 


0.25 


0.40 0.15 0.55 0.35 1.7e8 


2.4e8 


1.5e8 


1.9e9 


Bar-Unstable ~ L* Disk: At Peak of Inflow 


barex(pk) 


10 a ,50 


0.25 


0.48 0.20 0.13 0.23 l.le9 


2.5e8 


5.2e8 


3.8e9 


Bar-Unstable ~ L» Disk: After Bar Relaxation 


barex(re) 


10", 50 


0.25 


0.04 0.11 0.16 0.18 1.2el0 


5.5e9 


5.6e9 


1.5el0 



Parameters describing representative examples of our galaxy-scale simulations of galaxy-galaxy mergers and unstable isolated disks. The top row gives the 
range spanned in each parameter across our suite of simulations. Subsequent rows pertain to specific examples at chosen times in their evolution. (1) 
Simulation name/ID. (2) Minimum smoothing length (in pc). (3) Equation of state parameter (Figure 1). (4) Gas fraction = M gas / (M gas +M» t ji s i c ) of the 
disky/cold component (inside the given radius). (5) Scale height of the disky/cold component within the given radius (median \z\/R; the plane of the disk is 
defined by the total angular momentum vector). (6) Bulge-to-total mass ratio inside a given radius (bulge here includes all spherical components: stellar 
bulge, dark matter halo, and black hole). (7) Maximum surface density of the disky/cold component (gas plus stars) - for an exponential profile, the surface 
density is nearly constant at radii below the disk scale length. Otherwise, averaged over a couple times our minimum smoothing length. (8) Average surface 
density inside 300 pc for the bulge component. (9) Total mass enclosed inside a given radius. All simulations include black holes, but these are dynamically 
unimportant on these scales. 

" For each of the two simulations here, we have also run three ultra high-resolution simulations which also act as a moderate-resolution intermediate-scale 
simulations (~ 10 7 particles). They are identical in initial conditions to the standard merger and isolated disk run here, with initial gas fraction equal to, 
one-half, and one-quarter that shown (six simulations in total). 



Freeman 2004; Mandelbaum et al. 2006, and references therein). 
The galaxies have total masses = V v ] r /(WGH[z]) for an initial 
redshift z, with the baryonic disk having a mass fraction ma (typi- 
cally ma ~ 0.041) relative to the total mass. The system has an ini- 
tial bulge-to-total baryonic mass ratio B/T, and the disk has initial 
gas fraction / gas . The dark matter halos are assigned a concentration 
parameter scaled as in Robertson et al. (2006b) for the galaxy mass 
and redshift following Bullock et al. (2001). Disk scale lengths are 
set in accordance with the above scaling laws. 

In previous papers (referenced above), a large suite of these 
simulations have been presented, with several hundred simula- 
tions of varying equation of state, numerical resolution, merger 
orbital parameters, structural properties (e.g. profile shapes, initial 
bulge-to-disk ratios, and scale lengths), initial gas fractions, and 
halo concentrations. In this suite, galaxies have baryonic masses 
~ 10 s — 10 13 Mq and gas fractions / gas =0—1; mergers spanning 
mass ratios from 1:1 to 1:20, and isolated disks have Toomre Q 
parameters from 0.1 — 10. 

In this work, we focus on galaxies with baryonic masses 
~ 10 u Mq. Based on the survey above, we select a representa- 
tive simulation of a gas-rich major merger and one of an isolated 
disk, to provide the basis for our subsequent re-simulations. Some 
of the salient parameters of these simulations are given in Table 1 . 
The merger is equal-mass, with 10 11 Mq galaxies that have gas 
fractions of ~ 40% at the time of the merger/coalescence. The or- 
bit is a moderately tilted prograde, parabolic case (orbit e in Cox 
et al. 2006). Together this makes for a fairly violent, gas rich major 
merger, representative of many of our other gas-rich major merger 



simulations at both low and high redshifts. The isolated system is a 
10 11 Mq disk with / gas = 0.4, B/T = 0.3, scale length h = 3.2kpc, 
and mj — 0.041; it has a Toomre Q of order one and develops a 
moderate bar (amplitude ~ 15%), but the gas encounters an inner 
Linblad resonance at ~l-2kpc. 

Small variations in the orbits or the structural properties of the 
galaxies will change the details of the tidal and bar features on large 
scales. However, the precise details of these large-scale simulations 
are not important for our study of the dynamics on small scales (see 
Appendix A). Rather, the small-scale dynamics depends on global 
parameters such as the gas mass channeled into the central region, 
relative to the pre-existing bulge, disk, and black hole mass. In these 
respects, we have chosen the simulations summarized in Table 1 to 
be representative of a broad class of gas-rich systems. 

Our galaxy-scale simulations have spatial resolution - grav- 
itational softening length and minimum adaptive SPH smoothing 
length - of 50pc. In the suite described above, the resolution scales 
with galaxy mass and is ~ 50— 100 pc for M„ ~ 10" Mq systems, 
but in a subset of higher-resolution cases is as small as 20pc. In 
Hopkins et al. (2008b) and Hopkins et al. (2009a) we have demon- 
strated that this resolution is sufficient to properly resolve not only 
the mass fractions but also the spatial extent of the "starburst" 
formed from gas which loses angular momentum in a merger or 
via a strong bar instability. However, to assess how much of this 
gas can ultimately fuel a central BH requires that we determine the 
dynamics on even smaller spatial scales. 
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2.3 Intermediate Scales: Re-Simulating from 1 kpc to 10 pc 

In order to follow the behavior of gas inflow on smaller scales, we 
re-simulate the central regions of interest at higher resolution, in a 
series of progressively smaller-scale runs. We begin by selecting a 
number of representative outputs from the galaxy-scale simulations 
described above, near the peak of activity. We select several snap- 
shots in the gas-rich merger at key epochs: early close encounters 
of the two galaxies, just at nuclear coalescence (which is the peak 
of star formation in the nuclear region), and at the "end" of the 
merger (roughly ~ 10 8 yr after the final coalescence). We also se- 
lect snapshots typical of isolated, moderately bar-unstable systems, 
at times where a bar and some inflow has developed; for compar- 
ison, we also consider a fully stable (pre-bar) galaxy disk. In each 
case, we focus on the central 0.1 — 1 kpc region, which includes the 
majority of the gas that has been driven in from larger scales. Some 
of the representative properties of these snapshots, at these scales 
and times, are outlined in Table 1. 

Our approach to re-simulating the nuclear region is to use 
the larger-scale simulations to motivate the initial conditions of a 
smaller scale calculation (a "zoom-in" or "re-simulation"). We do 
so by de-composing the potential, density, and velocity distribu- 
tions of the gas, stars, and dark matter at a given time in the larger- 
scale simulation using the basis expansion proposed in Hernquist & 
Ostriker (1992). This allows us to not only re-construct a smoothed 
density profile, but also to include the asymmetric structures from 
the larger-scale simulation (if desired) and to define where the po- 
tential is noise-dominated. 1 From these stellar, gas, and dark matter 
distributions, we re-populate the gas and stars in the central regions 
(the scale we wish to re-simulate; generally out to an outer radius of 
~ 1 — 2 kpc) and use this as the initial condition for a new simula- 
tion that we run for several local dynamical times. To be conserva- 
tive, we typically initialize only a small amount of gas in the inner 
parts of the re-simulation 2 , since the larger-scale simulation from 
which the initial condition is drawn has little information about the 
gas properties on small scales; in Appendix A we show that the 
subsequent dynamics does not depend significantly on these details 
of the initial conditions. 

We have carried out a total of ~ 50 simulations at these inter- 
mediate scales, which together span a wide range in the key param- 
eters: the equation of state of the gas and the relative mass fraction 
in a pre-existing bulge, gas disk, and stellar disk. Table 2 summa- 
rizes the key properties of physical importance in several of these 
simulations (some numerical studies and surveys of initial condi- 
tion, which turn out to have little effect on the key results, are not 
listed). Because our general approach is to systematically survey 
the initial conditions, we do not identify every simulation in Table 2 
with an exact snapshot from Table 1 ; rather, they should be consid- 
ered a systematic parameter survey of possible intermediate- scale 
conditions, motivated by the typical range of sub-kpc conditions 
seen in our galaxy-scale simulations. Dark matter is present, but is 



As one continues the expansion to include arbitrarily small-scale modes, 
the best-fit mode amplitude will eventually yield an amplitude consistent 
with the shot noise in the simulation, roughly at the scale of the median 
inter-particle spacing; we discard higher mode numbers as they are particle- 
noise dominated. 

2 To avoid numerical effects from a step-function cutoff in the mass profile 
at small radii, we typically truncate the gas mass profile with a E gas oc 
(R/Ro) 2 power law inside of a radius So ~ 3 — 5 e in the parent simulation 
(e is the minimum smoothing length). Gas within this radius (< 1% of the 
re-simulated mass) is initialized with circular orbits. 



dynamically irrelevant at these scales. We have also varied the mass 
of the central black hole, but at these scales it is still dynamically 
unimportant. Although the initial conditions for our calculations 
are drawn from galaxy-scale simulations, the dynamics on small- 
scales depends primarily on a few key properties of the simulation 
(/gas and B/T), and is thus insensitive to many of the details of the 
galaxy on larger scales. 

Our intermediate scale simulations typically involve ~ 
10 6 particles, with a force resolution of a few pc and a particle mass 
of « 10 4 Mq. The duration of the simulation is ~ 10 7 - 10 8 yr - 
this is many dynamical times at small radii, but small compared to 
the dynamical time at larger radii. 3 These re-simulations can thus 
be thought of as a probe of the instantaneous behavior of the gas at 
small radii given the inflow conditions set at larger radii. 

We discuss a number of numerical tests and variations about 
this basic methodology in Appendix A. Specifically, we show that 
our results are not sensitive to including properties of the larger- 
scale galaxy in which the simulation should be embedded, such 
as the ~ kpc-scale tidal potential. They are also not sensitive to 
whether we initialize an axisymmetric gas/potential distribution, or 
whether the initial condition includes the non-axisymmetric modes 
present in the larger-scale simulation. The reason is that instabili- 
ties due to self-gravity can grow exponentially from shot-noise in 
the simulation even given an initially axisymmetric structure. Thus 
the presence of initial asymmetries on ~ 100 pc scales does not 
have a significant effect on the resulting transport of gas to smaller 
radii; the transport is determined by the presence (or absence) of 
internal instabilities in the gas on small scales. Finally, because the 
"initial" densities on small scales are intentionally initialized to be 
low relative to their values later in the simulation, our results are 
not particularly sensitive to how we initialize gas at small radii in 
the re-simulation; this is important because there is no reliable in- 
formation about these small-scales in our larger-scale simulation. 

To check that our re-simulation approach has not introduced 
any artificial behavior, we have run a small number of ultra-high 
resolution galaxy scale simulations, the inner properties of which 
can be compared to the intermediate scale re-simulations summa- 
rized here. For these very high resolution calculations there is con- 
tinuous inflow from large scales, so they can be self-consistently 
evolved for many dynamical times. The expense of these cal- 
culations, however, limits the survey of initial conditions possi- 
ble. We have 6 such simulations: three mergers, identical in mass 
and geometry to our canonical case in Table 1, but with initial 
/gas = 0.2, 0.4, 0.8, and ~ 10 7 particles, which gives SPH smooth- 
ing lengths of ~ lOpc. While not quite as high-resolution as our 
re-simulation runs, these provide an important check on the results 
of the latter and are run self-consistently for 4 x 10 9 yr. We will 
show results from these ultra-high resolution simulations at several 
points; we find that they are quite similar to our re-simulation runs, 
thus supporting the methodology used for most of our calculations. 

The very high resolution merger calculations also allow us to 
follow the binary BH pair to much smaller separation (our assump- 
tions lead to rapid merger below « lOpc in these calculations). We 
confirm, in these cases, that the BH-BH merger precedes most of 
the gas inflows at > 10 pc, so that the assumption that the binary 
has merged is probably reasonable for our re-simulation calcula- 



3 To ensure there are no later-time phenomena of interest, and to study the 
relaxed structure of the stellar remnant produced by each re-simulation, we 
evolve most for > 2 X 10 s years, by which time all the gas is exhausted. We 
find that there is no qualitatively new behavior at these later times. 
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Table 2. Intermediate-Scale Resimulations (~ 10 — 1000 pc) 



Simulation 


e 


#eos 


/gas 


h/R 


B/T(<R) 


S d (0) 


S b (0) 


M t0 ,(<R) [M ] 


Name 


[pc] 




(500 pc) 


(500 pc) 


lOOpc 


300 pc 


[M kpc- 2 ] 


lOOpc 


300 pc 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 




(7) 


(8) 




(9) 


If9bS" 


1.0 


b ,0.175,0.25 


0.90 


0.30 


0.5 


0.15 


l.OelO 


l.lelO 


5.4e8 


2.9e9 


If9bSthin 


1.0 


0.125,0.25 


0.90 


0.08,0.16 


0.4 


0.2 


l.OelO 


l.lelO 


5.4e8 


2.9e9 


If9b5res 


0.3,1,3,10 


0.125 


0.90 


0.30 


0.5 


0.15 


l.OelO 


l.lelO 


5.4e8 


2.9e9 


If9b5q 


1.0 


0,0 C ,0. 125,0.25,0.5,1 


0.90 


0.30 


0.5 


0.15 


l.OelO 


l.lelO 


5.4e8 


2.9e9 


Ilowresq 


3.0 


0,0.25,1 


0.95 


0.27 


0.0 


0.0 


6.0el0 


0.0 


1.2e9 


4.9e9 


Iflbllate 


1.0 


0.125 


0.096 


0.25 


0.06 


0.03 


l.Oell 


l.lelO 


3.4e9 


2.2el0 


IflbOlate 


1.0 


0.25 


0.091 


0.28 


0.002 


0.003 


6.0el0 


5.0e8 


1.3e9 


5.3e9 


IflbOlatel/ 


2.0 


0.25 


0.091 


0.28 


0.005 


0.01 


6.0el0 


5.0e8 


7.4e9 


8.5e9 


If3b3mid 


1.0 


0.125 


0.34 


0.25 


0.3 


0.15 


3.1el0 


2.0el0 


1.3e9 


7.2e9 


If3b3midRg< 


1.0 


0.125 


0.45,0.20,0.05 


0.3,0.5 


0.26 


0.12 


3.6el0 


2.0el0 


1.5e9 


9.2e9 


Iflb3Lmid 


1.0 


0.25 


0.10 


0.30 


0.07 


0.10 


6.4el0 


1.6e9 


1.4e9 


5.7e9 


Iflb3LmidL 1 ' 


2.0 


0.25 


0.10 


0.30 


0.15 


0.25 


6.4el0 


1.6e9 


8.4e9 


l.lelO 


If5b4mbul 


1.0 


0.25 


0.50 


0.24 


0.40 


0.55 


7.4e8 


1.9e8 


().3e8 


1.3e8 


If5b8mbul 


1.0 


0.25 


0.50 


0.24 


0.80 


0.90 


7.4e8 


3.1e9 


1.2e8 


5.3e8 


If9bllowm 


1.0 


0.25 


0.90 


0.15 


0.03 


0.06 


6.4e9 


1.2e7 


1.3e8 


5.4e8 


If3b9dsk 


1.0 


0.25 


0.32 


0.22 


0.95 


0.80 


1.6e9 


l.lell 


2.1e9 


6.2e9 


If3b9dskL rf 


2.0 


0.25 


0.32 


0.22 


0.71 


0.62 


1.6e9 


l.lell 


8.0e9 


9.3e9 


IfXb2gas 


1.0 


0.25 


0.16,0.32,0.50 


0.26 


0.16 


0.08 


3.6el0 


l.OelO 


1.3e9 


7.7e9 


Inf28b2 


1.0 


0.20 


0.20,0.80 


0.25 


0.51 


0.28 


6.3e9 


3.0e8 


3.6e8 


1.7e9 


Inf28b4 


1.0 


0.20 


0.20,0.80 


0.25 


0.67 


0.43 


3.3e9 


3.0e8 


2.7e8 


l.le9 


Inf28b6 


1.0 


0.20 


0.20,0.80 


0.25 


0.78 


0.56 


3.3e9 


6.0e8 


4.0e8 


1.4e9 


Inf28b8 


1.0 


0.20 


0.20.0.80 


0.25 


0.92 


0.81 


1.0e9 


6.0e8 


3.4e8 


9.5e8 


Inf2b9 


1.0 


0.20 


0.20 


0.25 


0.96 


0.89 


5.2e8 


6.0e8 


3.2e8 


8.6e8 


Inf28b2b f 


0.3 


0.125 


0.20,0.80 


0.25 


0.38 


0.21 


5.5e9 


5.0e7 


2.3e8 


Ue9 


Inf8b2h/ 


0.3 


0.20 


0.80 


0.25 


0.17 


0.11 


5.5e9 


2.5e7 


2.2e8 


l.le9 


Inf28b9b f 


0.3 


0.125 


0.20,0.80 


0.25 


0.81 


0.75 


5.5e9 


5.0el0 


2.1e9 


l.lelO 



Parameters describing our re-simulations of the 0.01 — 1 kpc regions from galaxy scale simulations. Parameters separated by commas denote 
simulations with otherwise identical initial conditions, re-run with the specified parameter varied. (1) Simulation name/ID. (2) Minimum 
smoothing length (in pc). (3) Equation of state parameter (Figure 1). (4) Initial gas fraction of the disky/cold component (inside the given radius). 
(5) Initial scale height of the disk component (inside the given radius). (6) Initial bulge-to-total mass ratio inside a given radius (again, bulge 
refers to all spherical components). (7) Initial maximum surface density of the disky/cold component (gas plus stars). (8) Initial average surface 
density inside 300 pc for the bulge component (9) Initial total mass enclosed inside a given radius. All simulations include BHs and dark matter, 
but these are dynamically unimportant on these scales. 

" A series of 7 runs testing different means of constructing initial conditions, described in Appendix A. 
* Isothermal equation of state, but with a large c s = 50kms~' cooling "floor." 
c Cooling allowed down to 100 AT, i.e. c s = 1 kms" 1 . 

d Somewhat larger-scale simulation (between "galaxy scale" and standard "intermediate scale"). Instead of B/T(< R) and M lot (< R) being 

evaluated at 100 pc and 300pc, they are here evaluated at 500 pc and lkpc, respectively. 
e Series where the gas disk profile is allowed to vary independent of the stellar disk profile. The gas has exponential, power-law, and truncated 

power-law profiles, with varying concentrations with respect to the disk (for example including an extended gas "reservoir" at a distance ~ 2 

times the regular nuclear stellar disk length, with surface density profile £ <x R~ '). 
f Very high-resolution simulations which also act as a moderate-resolution nuclear-scale simulations (2 X 10 7 particles; gas particle mass 

50()Mq). A series of 6 galaxy-scale runs with very high (~ lOpc) resolution, used as moderate-resolution intermediate-scale simulations, are 

also described in the text. 



tions. Moreover, the gas mass at ~ 10 pc is large (~ Mbh) in the 
merger simulations. Thus if gas-rich reservoirs indeed drive rapid 
BH-BH coalescence, the rapid merger of the two BHs should be 
a reasonable assumption on all of the scales that we simulate (see 
e.g. Escala et al. 2004; Perets et al. 2007; Mayer et al. 2007; Perets 
& Alexander 2008; Dotti et al. 2009; Cuadra et al. 2009). 

2.4 Nuclear Scales: From 10 pc to 0.1 pc 

The characteristic initial scale-lengths of the nuclear disks in our 
intermediate scale calculations are ~ 0.2 — 0.5 kpc. As we discuss 
in §3, if the gas fraction is sufficiently large, instabilities quickly 
develop that transport material down to ~ 1 — lOpc, near the reso- 
lution limit of our intermediate scale calculations. Material begins 



to pile up at these radii because the BH mass dominates the po- 
tential and the efficiency of large-scale modes decreases at small 
radii. In order to understand the dynamics on yet smaller scales, we 
therefore repeat our "re-simulation" methodology once more. The 
approach is identical to that described above, but this time using 
the intermediate-scale simulations with resolution of < lOpc as our 
"parent simulation" from which to motivate the initial conditions. 

We again carried out ~ 50 such simulations, typically with 
~ 10 6 particles, and a force/spatial resolution of ~ 0.1 pc (particle 
mass « IOOMq). The properties of these simulations are summa- 
rized in Table 3. The simulations are evolved for ~ 10 5 — 10 6 yr; 
this is large compared to the dynamical time at the smallest radii 
~ 0. 1 pc, but very small relative to the dynamical time of the larger- 
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Table 3. Nuclear-Scale Resimulations (~ 0.1 - 100 pc) 



Simulation 


e 


#eos 


/gas 


h/R 


M B h 


M b/a 


(<ff) 


S d (0) 


M d (<R) [M ] 


Name 


[pc] 




(100 pc) 


(100 pc) 


[Mq] 


10 pc 


50 pc 


[M kpc- 2 ] 


1 pc 


lOpc 


50 pc 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 




(9) 




Nf8hlc0 


0.3 


0.0,0.25 


0.75 


0.16,0.26 


2.9e7 


1.8e5 


3.0e5 


7.5el0 


1.0e5 


1.5e6 


1.0e8 


Nf8hlc0hol fl 


0.3 


0.25 


0.75 


0.24 


2.9e7 


1.8e5 


3.0e5 


7.4el0 


1.0e3 


1.2e6 


1.0e8 


Nf5h2cl 


0.1 


0.25 


0.50 


0.24 


3.0e8 


1.2e7 


8.0e7 


1.2ell 


3.6e5 


2.4e7 


1.6e8 


Nf28hlcl 


0.1 


0.25 


0.19,0.75 


0.27 


3.0e7 


1.3e7 


2.8e7 


7.5el0 


3.6e5 


2.4e7 


1.6e8 


Nf7hl2c0dsk 


0.1 


0.25 


0.70 


0.23 


2.9e7,3.0e8 


0.0 


0.0 


4.5e9 


1.0e4 


0.9e6 


6.0e6 


NfShlcO 


0.1 


0.25 


0.48 


0.25 


2.9e7 


0.0 


().() 


1.2ell 


3.0e5 


2.4e7 


1.6e8 


Nf5hlc2 


0.1 


0.25 


0.48 


0.25 


2.9e7 


3.0e7 


9.0e7 


1.2ell 


3.0e5 


2.4e7 


1.6e8 


Nf8hlc0thin 


0.1 


0.125 


0.75 


0.16,0.27 


3.0e7 


1.8e5 


0.3e7 


7.6el0 


2.0e5 


1.5e7 


1.0e8 


Nf5hlclthin2 


0.1 


0.125,0.25 


0.50 


0.07,0.14 


3.0e7 


3.0e5 


1.2e7 


1.2ell 


2.0e5 


2.4e7 


1.6e8 


Nf8hlclqs 


0.1 


0,0.125,0.25,1 


0.75 


0.28 


3.0e7 


3.0e6 


1.4e7 


7.3el0 


2.0e5 


1.7e7 


1.2e8 


Nf8h2cl 


0.1 


0.25 


0.75 


0.2,0.33 


3.0e8 


3.6e6 


1.4e7 


7.6el0 


2.0e5 


1.9e7 


1.2e8 


Nflhlcllow 


0.1 


0.25 


0.08 


0.25 


3.0e7 


3.5e6 


1.4e7 


1.7ell 


4.7e5 


3.7e7 


2.5e8 


NBhlclmid 


0.1 


0.20 


0.26 


0.23 


3.0e7 


3.5e6 


1.4e7 


2.1ell 


6.0e5 


4.6e7 


3.0e8 


Nf6hl2c2dsk 


0.1 


0.25 


0.57 


0.22 


2.9e7,3.0e8 


7.2e6 


2.6e7 


4.4e9 


l.le5 


8.1e6 


3.4e7 


Nf8hlc3dskM 


0.1 


0.25 


0.75 


0.30 


2.9e7 


1.6e7 


7.1e7 


7.4el0 


3.5e5 


3.0e7 


1.7e8 


NfShlcldens 


0.1 


0.25 


0.75 


0.25 


3.0e7 


3.6e6 


1.4e7 


3.8ell 


Ue6 


8.1e7 


5.3e8 


NfShlclICs* 


0.1 


0.25 


0.75 


0.28 


3.0e7 


3.1e6 


1.4e7 


7.3el0 


2.0e5 


1.7e7 


1.2e8 


Nf8hlclthin 


0.1 


0.125,0.18 


0.75 


0.08,0.17 


3.0e7 


3.1e6 


1.4e7 


7.3el0 


2.0e5 


1.7e7 


1.2e8 


Nf2h2b2 


0.1 


0.20 


0.20 


0.25 


3.0e7 


6.5e6 


2.0e7 


3.0ell 


8.1e4 


7.0e7 


4.2e8 


Nf8h2b2 


0.1 


0.20 


0.80 


0.25 


3.0e7 


1.3e7 


4.0e7 


1.5ell 


4.7e4 


3.5e7 


2.1e8 


Nf2h2b4 


0.1 


0.20 


0.20 


0.25 


3.0e7 


6.5e6 


2.1e7 


1.5ell 


4.3e4 


3.5e7 


2.1e8 


Nf8h2b4 


0.1 


0.20 


0.80 


0.25 


3.0e7 


1.3e7 


4.0e7 


7.7el0 


2.4e4 


1.7e7 


l.le8 


JNf2h2b5 


0.1 


0.20 


0.20 


0.25 


3.0e7 


6.5e6 


2.1e7 


7.7el0 


2.1e4 


1.8e7 


l.le8 


Nf28h2b6 


0.1 


0.20 


0.20,0.80 


0.25 


3.0e7 


l.le7 


3.7e7 


3.8el0 


l.le4 


8.8e6 


5.3e7 


Nf8h2b8 


0.1 


0.20 


0.80 


0.25 


3.0e7 


1.3e7 


4.1e7 


1.5el0 


4.7e3 


3.5e6 


2.1e7 


Nf28h2b9 


0.1 


0.20 


0.20,0.80 


0.25 


3.0e7 


l.le7 


3.8e7 


9.6e9 


2.7e3 


2.2e6 


1.3e7 


Nf8h2blh 


0.015 


0.20 


0.80 


0.25 


3.0e7 


6.4e6 


2.0e7 


1.5ell 


5.2e4 


3.5e7 


2.1e8 


Nf8h2b3L 


0.1 


0.20 


0.80 


0.25 


3.0e7 


2.3e7 


l.le8 


1.5ell 


9.3e3 


3.9e7 


4.9e8 


Nf8h2b4q 


0.1 


0,0.02,0.06,0.12, 
0.25,0.35,0.5,0.7,1 


0.80 


0.25 


3.0e7 


1.3e7 


4.0e7 


7.7el0 


2.4e4 


1.7e7 


l.le8 



Parameters describing our nuclear-scale re-simulations of the sub-100 pc regions from intermediate-scale simulations. (1) Simulation name/ID. (2) 
Minimum smoothing length (in pc). (3) Equation of state parameter (Figure 1). (4) Initial gas fraction of the disky/cold component. (5) Initial scale 
height of the disky component. (6) Black hole mass (Mq). (7) Initial bulge or nuclear stellar cluster mass, inside the given radius. (8) Initial 
maximum surface density of the disky/cold component (gas plus stars). (9) Initial mass of the disky component (gas plus stars) inside a given radius 
(does not include the BH mass or, if significant, nuclear star cluster/bulge mass). Dark matter is insignificant on these scales. 

" Central "hole" is extended in disk out to 5 pc. 

* Simulations with no central deficit of matter; the initial density from the larger-scale simulation is extrapolated in to r — > 0. Also expanded into series 
of initial conditions, as described in Appendix A. 



scale simulations from which the initial conditions are drawn. The 
characteristic spatial scale of the re-simulated material is initially 
~ 10 — 30pc. As described in Appendix A, we carried out a num- 
ber of numerical tests of the robustness of these simulations. 

At radii ~ pc, the parameters that determine the dynamics are 
largely the equation of state of the gas, the mass of the BH, the 
mass of the nuclear disk formed by the inflow from larger scales, 
and the gas fraction of that nuclear disk. Since the BH dominates 
the spherical component of the potential at these radii, the "bulge" 
mass at these radii is only of secondary importance; we include it 
but find that it makes little difference. 

As in §2.3, we have checked the results of these "re- 
simulations" by carrying out a small subset of ultra-high resolu- 
tion runs. These extend from ~ 0.3 — 1000 pc and follow inflow 
from larger scales deep into the potential of the BH; because they 
resolve larger spatial scales than our typical "nuclear scale" simula- 
tion, these can be run self-consistently for 2 x 10 s yr. Specifically, 
we have five such high-resolution intermediate scale simulations 



(see Table 2), three with initial / gas = 0.8 (a low, intermediate, and 
high B/T case), and two with / gas = 0.2 (low and high B/T). They 
have ~ 10 7 particles and gravitational softening lengths of ~ 0.3 pc. 
We show the results from these runs explicitly at several points; 
we find that they are completely consistent with our survey of re- 
simulations, which cover a larger parameter space of galaxy/BH 
properties, but are more limited in dynamic range. 

It is important to note up-front that our simplified treatment of 
the ISM physics becomes particularly suspect on nuclear scales ~ 
pc. At these radii, our assumption that we can average over the dy- 
namics of stellar winds, supernovae, HII regions, etc. and define an 
effective ISM equation of state may break down. Nonetheless, we 
believe that the efficient angular momentum transport found here is 
likely generic, so long as some of the gas is prevented from forming 
stars and the gas fraction is sufficiently high that instabilities gen- 
erated by self-gravity are initiated. The fact that the main sequence 
lifetime of a massive star is longer than the local dynamical time 
on small scales probably increases the efficacy of stellar feedback 
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and decreases the fraction of the gas turned into stars per dynamical 
time (Murray et al. 2009). 

At scales <C 0. 1 pc, the potential is fully Keplerian and viscous 
heating is sufficient to stabilize the disk against its own self-gravity 
(i- e -> 2^1) (Goodman 2003). At these radii, the system begins to 
approach a traditional accretion disk. Given the cessation of star 
formation and the deep potential well of the BH, we assume that 
the inflow rate at ~ 0. 1 pc is a reasonable proxy for the true accre- 
tion rate onto the BH. Because our simulations are not well-suited 
to describe the physics of the disk on scales < 0. 1 pc, we do not 
perform a further "zoom in." 

3 OVERVIEW: FROM KPC TO SUB-PC SCALES 

Using the numerical simulations described in §2, we now describe 
how gas is transported from 3>kpc scales to <Cpc scales. Initially, 
our discussion is somewhat qualitative; we focus on emphasizing 
the key physics at play and our key results. In § 5 we discuss the 
relevant stability criteria more quantitatively, and outline some spe- 
cific criteria necessary for "interesting" gas inflow. 

Figure 2 shows an illustrative example of the results of our 
re-simulations on various scales. We plot gas surface density maps, 
with color encoding the gas effective sound speed, from scales of 
~ 100 kpc to < 1 pc. The simulation in this case is a fairly gas-rich 
major merger (/ gas ~ 30 — 40% at the time of the final coalescence) 
of two 5 x 10 10 Mq baryonic mass galaxies. The smaller scale re- 
simulations were carried out just after the coalescence of the two 
nuclei, which is near the peak of star formation activity, but when 
the system is still quite gas rich. The initial systems had pre-existing 
bulges of ~ 1/3 the disk mass and BHs initialized on the Mbh — <J 
relation (~ 10 7 Mq). Each panel is rotated so that the view is close 
to "face-on" with respect to the total angular momentum vector 
of the gas plotted in the image. Viewed edge-on, much of the gas 
forms a modestly thin (H/R < 0.3) disky distribution at all radii. 

3.1 Large Scales: Mergers and Bars 

From ~ 100 kpc to > lOOpc, the gas flows are well-resolved by 
our galaxy-scale simulation. In mergers, the final collision of the 
two galaxies yields strong torques that efficiently cause most of the 
gas to flow to the center on a timescale approximately equal to a 
few dynamical times. This process has been described in detail in 
e.g. Hopkins et al. (2009b), and references therein, but for com- 
pleteness we briefly summarize the important physics. The sec- 
ondary/merging galaxy does not directly torque the gas. Rather, 
the torques on the gas are dominated by local torques from stars 
originally in the same disk as the gas. The merger induces triaxial 
"sloshes" and bar-like structures in the stars, i.e., non-axisymmetric 
modes, supported by radial and/or random orbits. These are also 
induced in the gas. However, because the gas is dissipational, the 
gaseous modes slightly lead those in the stars. The stellar distur- 
bance, being physically close to, trailing, and in near-resonance 
with the gas, produces a strong torque that removes angular mo- 
mentum from the gas and easily dominates the total torque (torques 
directly from the secondary galaxy, from the primary halo, or hy- 
drodynamic torques from shocks or internal clump collisions, are 
all < 10% effects; see e.g. Barnes & Hernquist 1996; Barnes 1998; 
Hopkins et al. 2009b). 

After the two galactic nuclei coalesce, the disturbances in the 
stellar component of the galaxy relax away in a number of crossing 
times. Until they relax, gas inflows continue. Moreover, the coales- 
cence of the nuclei completes much more rapidly, in a timescale 
close to a single crossing time, at least in a major merger. Thus, a 
significant fraction of the gas inflows can occur in the background 



of a rapidly relaxing stellar potential, in the wake of the nuclear 
coalescence. This is the stage illustrated in Figure 2. 

The gas that loses angular momentum flows in to radii ~ 
0.5 — 1 kpc (for an ~ L» system), where it participates in a nuclear 
starburst and builds a dense central stellar mass concentration, criti- 
cal for establishing the structural properties and size of the remnant 
spheroid. At these scales, the system is often gas-dominated for a 
short period of time owing to these inflows (provided the merger 
is sufficiently gas-rich). However, as the gas forms stars, the cen- 
tral region will quickly become more stellar-dominated; because 
these stars form out of the gaseous disk, in the relaxing potential 
— they are not themselves violently relaxed. This is important for 
the subsequent evolution of the system because of the presence of 
disk instabilities that would be suppressed by a larger dispersion- 
supported (spherical) component in the very central region. 

The general scenario summarized here can be applied not just 
to major mergers, but also to minor mergers, fly-by encounters, and 
even sufficiently bar-unstable stellar disks. The details will be dif- 
ferent, but the qualitative steps above, and the exchange of angular 
momentum between gas and stars, is robust, ultimately leading to 
inflow to sub-kpc scales. The subsequent evolution depends largely 
how much material is efficiently channeled to small radii (relative 
to the bulge and BH mass), not on how that material gets there. 

3.2 Morphology and Gas Transport From 1 kpc to 10 pc 

3.2.1 General Behavior 

The gas Mailing from large radii begins to "pile up" at radii 
~ 0. 1 — 1 kpc, rather than continuously flowing in to yet smaller 
radii. This is because the torques from the disturbances at large 
radii become less efficient at small radii. This happens for three 
reasons: (1) In the merger context, the stellar perturbation at small 
radii relaxes after coalescence, decreasing the efficiency of gas in- 
flow. (2) The rapid gas inflow implies that the system becomes in- 
creasingly gas-dominated at radii ~ 100 pc, even if the initial disk 
gas fraction is low, ~ 0. 1 . Because the primary angular momentum 
sink of the gas is the local stars, when the system becomes locally 
gas-dominated, angular momentum transfer is actually less efficient 
(see Hopkins et al. 2009b,f). (3) The gas can encounter the equiv- 
alent of an inner Linblad resonance. This is especially important 
for unstable gas bars, minor mergers, and disturbances induced by 
early passages. For the case of coalescence following major merg- 
ers, the disturbance is not a single mode, but a series of modes at all 
scales. As such, there is often no formal inner Linblad resonance or 
angular momentum barrier (each mode may have such a barrier, but 
these are spread over a wide range of scales; there is thus a means 
to overcome the barrier associated with any single mode). 

Figure 2 shows the outcome of gas pile-up in the central 
kpc using an intermediate- scale re-simulation (middle row). In this 
case, the intermediate scale simulation is a high-resolution re- 
simulation of the larger-scale gas distribution at a given epoch in 
a gas-rich major merger. The gas density reached from the larger- 
scale inflows is quite large - ~ 10 10 Mq worth of gas has formed 
a disky component with a scale length of 0.3 kpc and an average 
surface density of ~ 10 10 Mq kpc -2 . This is a large fraction of the 
galaxy mass - larger than the pre-existing bulge within these radii. 
The small-scale gas disk is therefore strongly self-gravitating. In- 
deed, we see from Figure 2 that it quickly develops unstable, non- 
axisymmetric modes. 4 This is essentially the "bars within bars" 

4 These modes develop almost identically even if we initialize the re- 
simulation to be perfectly smooth and remove all external tidal forces (see 
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Figure 2. Example of our multi-scale simulations used to follow gas flows from ~ 100 kpc to ~ 0.1 pc. Each row is a separate simulation, with the initial 
conditions for the intermediate and nuclear-scale simulations taken from the output of the larger-scale runs in the row above it. Each panel shows the projected 
gas density (intensity) and effective sound speed (color; blue is gas with an effective c s ~ 10kms~', through yellow at ~ 100 — 200kms~ '). Each image 
is rotated to project the gas density "face on" relative to its angular momentum vector. From top left to bottom right, panels zoom in to the nuclear region 
around the BH, with resolution spanning a factor ~ 10 6 in radius. Top: Large-scale gas-rich galaxy-galaxy major merger simulation, just after the coalescence 
of the two nuclei (run b3ex(co) in Table 1). The apparent second nucleus is actually a clump formed from gravitational instability. Middle: A higher-resolution 
re-simulation of the conditions in the central kpc (run If3b3midRg in Table 2). Despite the fact that the background potential is largely relaxed on these scales, 
the very large gas inflows lead to a strongly self-gravitating disk on ~ 0.5 kpc scales that develops a strong spiral instability, leading to efficient angular 
momentum transport to ~ lOpc. Again, some clumping appears (there is only one nucleus). Bottom: High resolution re-simulation of the central ~ 30 pc of 
the intermediate-scale simulation, with a resolution ~ 0. 1 pc (run Nf5hlc2 in Table 3). The potential is quasi-Keplerian, suppressing traditional bar and spiral 
instabilities, but the large inflows lead to a self-gravitating system that develops a standing eccentric disk mode (single-armed m = 1). The stellar and gaseous 
eccentric disks precess relative to one another on ~ 1 — lOpc scales and drive efficient inflows of ~ IOMq yr _ 1 into the central 0. 1 pc. 
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Figure 3. Images of the morphologies of the secondary instabilities that develop on intermediate scales (~ 100 pc), as a consequence of galactic-scale torques 
bringing gas into the central kpc. Each panel shows projected gas density (intensity) and local star formation rate (color, from blue through yellow) in the 
central kpc. We sort the simulations by the parameters that have the largest impact on the dynamics: the bulge-to-disk ratio (disk being both stellar and 
gaseous) within 300pc and gas fraction M gas / (M gas +M*. disk): from top - bulge-dominated systems with B/T > 0.8 to bottom - disk-dominated systems 
with B/T < 0.1, and from left — systems with / gas < 0.1 to right - systems with / gas > 0.8. Global instabilities become more prominent with decreasing B/T; 
local instabilities (clumping/fragmentation and/or more tightly wound spirals) become more prominent with increasing / gas . The simulations at fixed B/T 
and/or / gas are similar, but include variation in e.g. the initial mass profile shapes; this contributes to the dramatic diversity of morphologies on nuclear scales. 
Almost all systems that are not entirely bulge-dominated develop strong large-scale instabilities that efficiently transport gas to ~ lOpc. Observationally, these 
would be categorized as a mix of bars, spirals, nuclear rings, crossed rings, fragmentation/clump instabilities, and single or triple-armed spirals. 



© 0000 RAS, MNRAS 000, 000-000 



How Do Massive Black Holes Get Their Gas ? 



13 



scenario predicted by Shlosman et al. (1989), although the mor- 
phology of the system is clearly not a simple bar; this is an impor- 
tant point to which we return below. The strength of the modes that 
develop in this re-simulation depend on the fact that there is star 
formation in the gas - as the gas turns into stars in situ, those stars 
develop non-axisymmetric modes, and the two precess relative to 
one another. As in the galactic-scale torques discussed above, this 
produces particularly efficient angular momentum transport. These 
processes ultimately lead to a significant amount of gas flowing 
down to ~ lOpc. 

3.2.2 Diversity in Morphologies and Inflow Strengths 

For our fiducial parameterization of the ISM equation of state, the 
disk-to-bulge ratio and gas fraction on ~ 0.1 — 0.3 kpc scales have 
the largest influence on the dynamics and angular momentum trans- 
port in our intermediate-scale simulations (for discussion of the role 
of <7eos and sub-resolution physics, see Appendix B). Figure 3 il- 
lustrates this by showing the ~ lOOpc scale morphology of a repre- 
sentative subset of our intermediate scale re-simulations, each after 
a couple of local dynamical times of evolution (~ 10 7 — 10 8 yr). 
These are sorted by gas fraction f g and B/T. 5 

As Figure 3 shows, systems with very large B/T > 0.9 (top 
row) are globally stable, as expected analytically. In the extremely 
gas-rich, large B/T, case (top right), local Toomre-scale instabili- 
ties develop - if such small clumps were infinitely long-lived, their 
orbits would decay via dynamical friction and allow for some trans- 
port of gas to the center. However, because the clumps are dense, 
they quickly turn into stars - such a mechanism is largely a means 
to move stars to the galactic nucleus, not gas (and in fact leads to 
nuclear stellar clusters much larger than those observed; see Ap- 
pendix B for a more detailed discussion). Previous claims that such 
clump-sinking could efficiently fuel BH growth (see e.g. the dis- 
cussion in Wada et al. 2009; Kawakatu & Wada 2008, and refer- 
ences therein) have neglected star formation and have thus dra- 
matically over-estimated the inflow of gas via this process. Lo- 
cal clumping instabilities like these do, of course, occur, forming 
molecular clouds and star clusters. However, there are both theo- 
retical and observational arguments that suggest that they quickly 
dissolve on a few local dynamical times, probably via feedback 
from some combination of stellar winds, HII regions, and radiation 
pressure (Larson 1981; Blitz 1993; Krumholz et al. 2006; Allen 
et al. 2007; Bonnell et al. 2006; Elmegreen 2007). We also note 
that, under certain conditions, such local instabilities could instead 
produce small scale, tightly wound spiral waves instead of leading 
to fragmentation (e.g. Lodato & Rice 2004; Rice et al. 2005; Boley 
et al. 2006); however, this generally occurs when the cooling time 
is comparable to or larger than the dynamical time, which is not the 
case in these simulations (although it will also depend on the tur- 
bulent decay time which, if sufficiently large, can allow turbulence 
to suppress runaway clumping and star formation). 

Figure 3 demonstrates that the strength of global non- 
axisymmetric modes increases dramatically with decreasing B/T. 
In fact, as soon as the bulge and disk are comparable (second row 
from top), prominent modes appear. In the disk-dominated cases 
(bottom row), very large angular momentum transport occurs even 

Appendix A). It is thus not sensitive to the larger-scale environment; rather, 
the system is simply strongly globally unstable. 

5 The other parameters of the simulations are not all identical, representing 
the properties of the simulations from which they are selected, but the key 
qualitative behavior depends primarily on these two parameters. 



in cases without much total gas. Of course, at each B/T, increasing 
the gas content makes the system more vulnerable to local instabil- 
ities as well - usually making the overall inflow more clumpy and 
time-variable. In several cases, the large-scale mode (e.g. a spiral 
arm) becomes self-gravitating and globally fragments (not neces- 
sarily into smaller sub-units, but as a whole), leading to a major co- 
alescence - what is almost (dynamically speaking) a scaled-down 
merger in the central regions! 

Figure 3 also demonstrates an important point seen in Fig- 
ure 2. Although the instabilities seen in our simulations qualita- 
tively resemble the "bars within bars" idea of secondary instabili- 
ties once the gas density is sufficiently high, the morphologies vary 
widely, and are not restricted to traditional bars (although these cer- 
tainly do appear). At similar gas fraction and B/T, we find that the 
strength of angular momentum transfer is generally similar, but we 
also find that the visual morphology (and the precise modes im- 
portant for transport) can vary widely, depending on time and on 
the details of the gas, stellar disk, bulge, and halo profiles, and the 
precise equation of state. Thus, global quantities such as the mass 
profile and accretion rate are comparatively robust, but the obser- 
vational classification of these systems would vary widely. 6 Indeed, 
Figure 3 shows traditional bars and spirals, nuclear rings, crossed 
or barred rings, single or three-armed systems, flocculent disks, and 
clumpy, irregular morphologies. These all appear, with no obvious 
preference for one or another as a whole, in our simulations. 

This feature of our simulations may account for a number of 
observational results in the literature. For example, surveys of AGN 
have often found that although nuclear bars on these scales only ap- 
pear in some fraction of sources (not necessarily much larger than 
the fraction of non-AGN in which they appear), there are ubiquitous 
asymmetric gas structures of some sort, similar to those modeled 
here (see e.g. Martini & Pogge 1999; Garcfa-Burillo et al. 2005; 
Haan et al. 2009; Krips et al. 2007; Laine et al. 2002; Peletier et al. 
1999; Sakamoto et al. 1999). We discuss this further in § 7. 

3.3 Towards the Accretion Disk: Eccentric Disks at ~ Parsec 

From ~ 500 to ~ lOpc, our intermediate-scale simulations success- 
fully demonstrate efficient angular momentum transport via a wide 
range of gravitationally unstable modes. Near the smallest radii in 
these simulations, however, the systems encounter yet another an- 
gular momentum barrier. At that point, gas has reached the BH ra- 
dius of influence, i.e., the BH begins to contribute non-trivially to 
the potential, which becomes quasi-Keplerian. This halts further 
inflow because the disk is no longer strongly self-gravitating and 
so is less susceptible to global modes. In addition, the gas gen- 
erally encounters an inner Linblad resonance associated with the 
intermediate-scale bar. 

Indeed, it is widely appreciated that both "bars within bars" 
and the direct or induced torques due to perturbations from merg- 
ers, close passages, and large-scale bars do not produce efficient an- 
gular momentum transfer at radii ~ 0. 1 — lOpc (see e.g. the discus- 
sion in Athanassoula et al. 1983, 2005; Shlosman et al. 1989; Heller 
et al. 2001; Begelman & Shlosman 2009, and references therein). 



6 In detail, the angular momentum transport does depend on the structure of 
the unstable mode/perturbation. The precise dependence will be discussed 
in Paper II. At the dimensional level, the torques and inflow scale in the 
same manner independent of the detailed mode morphology, so long as the 
stellar torques are sufficiently strong to cause orbit crossing and shocks in 
the gas. The details of the specific modes driving such shocks amounts to 
numerical factors of ~a few in the torques. 
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Figure 4. Images of the instabilities that develop in our small-scale (~ 10 pc) nuclear re-simulations with 0.1 pc resolution. Each panel shows projected gas 
density (intensity) and local star formation rate (color, from blue through yellow) in the central kpc. Each of these simulations can be thought of as a simulation 
of the corresponding nuclear scales from Figure 3. The simulations extend into the BH radius of influence. The primary parameters of importance are the ratio 
of the gas mass to the BH mass (or BH plus bulge/star cluster mass, when the latter is present) and the gas fraction in the disky component: top to bottom 
is decreasing BH/stellar mass while the disk gas fraction increases from left to right. A strong m = 1 mode is generic for reasonable BH/stellar mass and 
gas fractions - this corresponds to an eccentric, globally precessing (non-winding) disk (or single-armed spiral), a mode that is special to the quasi-Keplerian 
potential. The resulting torques drive inflows of up to IOMq yr _ 1 at < 0. 1 pc scales (Figure 5), sufficient to fuel a luminous quasar. 



As a consequence this is often considered the most difficult-to- 
explain regime of gas inflow and angular momentum transport. 

We find, however, that efficient angular momentum transfer 
continues in our smaller-scale re-simulations for sufficiently gas- 
rich, disk-dominated systems. Figure 4 shows this with images 
spanning a range of f g and B/T; we will discuss this physics in 



more detail in §4 & 5. These results show that, so long as the gas 
density is sufficiently large (relative to the BH mass), the system 
develops a precessing eccentric disk (an m = 1 mode) that drives 
gas down to sub-pc scales ~ 0.1 pc. As before, stars rapidly form 
out of the disk, leading to a similar mode in both the stars and 
gas; these modes precess about the BH relative to one another with 
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Figure 5. Time-dependent inflow rate through various radii, for the example simulation shown in Figure 2. Main: Gas inflow rate into the central 300pc in 
the large-scale galaxy merger simulation, as a function of time since the beginning of the merger. The large spike follows the coalescence of the two nuclei. 
Dotted line shows the instantaneous inflow rates; solid line shows the inflow rate averaged over five local dynamical times. Top Right: Inflow rates into the 
central ~ 10 pc in the intermediate scale re-simulation for a small time interval just after coalescence. Top Left: Inflow rates into 0. 1 pc in our smallest-scale 
resimulation, initialized with conditions from the previous intermediate-scale simulation after it has evolved for several local dynamical times. The accretion 
rate reaches 1 — IOMq yr — , sufficient to fuel a luminous quasar. Note the episodic, highly variable nature of accretion/inflow at each scale. 



slightly different pattern speeds, leading to crossing orbits, dissi- 
pation of energy and angular momentum in the gas, and thus net 
inflow. This is, once again, an instability that depends primarily on 
the presence of sufficient gas at small radii in the first place. We 
find that this condition is met at some point(s) in time in all of our 
simulations with significant gas mass and instability at ~ lOOpc, 
i.e., in those simulations that meet the "bars within bars" criteria in 
§3.2 and Figure 3. 

4 INFLOW RATES AND GAS PROPERTIES 

Figure 5 shows the time-dependent inflow rates through several an- 
nuli, for the same set of multi-scale re-simulations shown in Fig- 
ure 2. These demonstrate and quantify the general scenario sum- 
marized above: on large scales, coalescence and the final stages of 
the merger drive a large quantity of gas into the central few hun- 
dred pc, with inflow rates ~ 100 — 300 Mq yr _1 . The total dura- 
tion of this phase is a few times 10 s yr. During this time, the gas 
accumulates at small radii ~ 100 pc; simulating the dynamics on 
this scale for ~ 10 7 yr, we find that secondary gravitational insta- 
bilities develop that drive further inflows into the central ~ lOpc, 
with inflow rates ~ 10 — 200 Mq yr~ 1 . Zooming in yet again dur- 
ing one epoch of significant inflow to < 10 pc, our smallest-scale 



simulations resolve the rapid formation of an eccentric nuclear disk 
around the central BH; the accretion rate into the central 0.1 pc, 
which is likely a reasonable proxy for the accretion rate onto the 
BH, reaches ~ 1 — lOMgyr -1 . This is sufficient to fuel a lumi- 
nous AGN at the Eddington rate. 

Figure 5 also demonstrates that the small-scale accretion rate 
can be highly time-variable. This is in part a consequence of the 
accretion of individual clumps/clouds (e.g., Fig. 4), but is also a 
consequence of the fact that gravitationally unstable perturbations 
rapidly grow, dissipate, and generate other structures; depending 
on, e.g., the state of precession of the stellar versus gaseous disk, 
the system can transition between inflow and outflow at a given 
radius. Even on the largest scales, the inflow is still highly vari- 
able, although is coherent over a time much longer than the dy- 
namical time because it is driven by the global torques involved in 
the merger. Because of the variability in M on different scales, we 
do not expect every merger (or isolated galaxy with a large-scale 
bar), at every time, to exhibit significant inflow from large scales 
all the way down to the BH. This is important - after all, a large 
fraction of observed mergers are not bright quasars. In addition, the 
large variation in the physical conditions at small radii for a rela- 
tively fixed set of conditions at larger radii demonstrates that great 
care must be taken when trying to correlate the galactic structure at 



© 0000 RAS, MNRAS 000, 000-000 



1 6 Hopkins and Quataert 




0.0 0.5 1.0 1.5 2.0 2.5 

t [Gyr] 



Figure 6. Time-dependent inflow rate through various radii, as in Figure 5, but for a simulation in which the large-scale inflow into the central kpc is weak 
(simulation barex in Table 1). Specifically, this is a gas-rich disk with B/T = 0.3, which is moderately bar-unstable. Compared to a major merger, the inflow 
to < 300 pc is much weaker, and the presence of an inner Linblad resonance at small radii leads to both outflow and inflow from the bar. Because of the much 
lower gas supply to the inner part of the galaxy, secondary instabilities do not efficiently develop that can transport gas further in. 



~ 0.1 — 1 kpc with the BH accretion rate in order to constrain the 
physics of AGN fueling. 

Figure 6 shows the same calculation of the inflow rates, but for 
re-simulations of an isolated disk galaxy with a moderate bar insta- 
bility (the case in Table 1). The system develops a fairly strong 
bar at about a kpc (representing a perturbation to the potential 
of (5<E>/$ ~ 0.15), but this is still much weaker than the major 
merger case shown in Fig 5; moreover, the pre-existing bulge with 
B/T = 0.3 means that there is an inner Linblad resonance at a 
few hundred pc. The net result is that the overall inflow is much 
weaker, and there can be outflow as well as inflow on small scales 
(the system overcomes the Lindblad resonance for short periods of 
time, leading to inflow followed by outflow as it re-equilibrates). 
Re-simulating smaller scales, the weak inflows produced by the 
bar lead to low gas mass (~ 1 — 10%) relative to the bulge mass 
on small scales, and global secondary instabilities do not develop. 
Inflow/outflow rates from the bar at 0.3 — 1 kpc are characteristi- 
cally ~ 0.5 — lOMQyr -1 , similar to observed local barred galax- 
ies (Quillen et al. 1995; Jogee et al. 2005; Haan et al. 2009); at 
~ lOpc they do not exceed lMQyr -1 , and at < O.lpc they are 
characteristically <C0.1 Mq yr~ 1 . At these low accretion rates our 
calculations are not reliable, but the key point is that in the low gas 
fraction limit, gravitational torques drive very little accretion from 
~ kpc to < 0.1 pc, even in the presence of a modest bar at ~ kpc. 



Figure 7 shows an example of the torques driving the inflows 
in our high accretion rate simulations. We show the radial torque 
profile at one time in one of our nuclear-scale simulations, broken 
down into gravitational torques from the stars and gas in the nuclear 
disk (versus the large-scale tidal field) and hydrodynamic torques 
(largely from pressure forces); note that negative torques drive in- 
flow. The average gravitational torque in an annulus between Ro 
and Ro + AR, from sources j, is defined by 



1 



- Mum J R[ 

1 



M g , s [AR] 



/%as(r) [r x a grav ,j] z drdcos6ld0 
^2 m dr t x V$ y (ri)], 



(2) 



where M gas [AW] = Yl m i i s me 8 as mass m me annulus; the sum 
over i in the second expression refers to the gas particles inside 
the annulus, while $ ; is the gravitational potential generated by 
the torquing particles of interest (gas, stars, etc.). 7 The radius r is 
defined with respect to the BH, and we focus on the z component of 



7 For small separations between two particles i and j, |rj — T\\ <e, we in- 
clude the GADGET force softening to be self-consistent (equivalent, inside 
a softening length, to the potential of a Plummer sphere), but this has little 
effect on our results. 
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Figure 7. Radial profile of the specific torque in a representative nuclear 
scale simulation (Nf8hlclqs in Table 3) during a stage of rapid BH ac- 
cretion. The torque per unit mass acting on the gas is averaged over narrow 
radial annuli at one time; negative torques decrease the gas angular momen- 
tum. The total gravitational torque (black line) is divided into the torque 
from stars in the same disk as the gas (purple), gas (red), and the contri- 
bution from the large-scale tidal field (orange). The hydrodynamic torques 
from pressure forces and artificial viscosity (green) are comparatively small. 
Inflow at all radii is dominated by gravitational torques, largely from nearby 
stars in the same m = 1 pattern as the gas. 



r, where the z axis is defined to be the angular momentum vector 
of the disk (torques in the other directions are negligible, so this is 
almost identical to plotting |r|). The torques from pressure forces 
are similarly defined by 



1 



TP A 



M gas [AR] 



R + AR 



Pgas(r) [r x ap] ; drdcos#d<^ 



1 



M m [AR] 



E m[r l x--rVP(r l )l 



(3) 



where the pressures and densities are determined in standard fash- 
ion from the SPH quantities. Note that the BH can and does move 
(with typical amplitude R ~ 0AR[M d (< R) ~ 0.1 M Bli ]) in re- 
sponse to the m — 1 modes on nuclear scales; but since we are in- 
terested in the inflow onto the BH itself, we evaluate these torques 
about its instantaneous position (rather than, say, the center-of-mass 
of the system). The results are qualitatively similar in Figure 7 if we 
fix the central position (albeit quantitatively altered within <C lOpc 
at a factor ~ 2 level), but are less directly relevant to interpreting 
inflow onto the BH. 

The qualitative behavior shown in Figure 7 is representa- 
tive of all of our simulations. At essentially all radii, gravitational 
torques dominate hydrodynamic torques. Moreover, the gravita- 
tional torques themselves are dominated by torques from stars, not 
the torques of the gas on itself; specifically, the stars that are impor- 
tant are in the same asymmetric perturbation as the gas. This is also 
the case on larger scales (> lOOpc), for both mergers and barred 
systems (see Barnes & Hernquist 1996; Hopkins et al. 2009b). 
Torques from the spherical component of the system (e.g. the halo 
and/or bulge stars) are negligible at all radii. Unsurprisingly, the 
torques from the large-scale tidal field, defined here as the torques 
from the 3> lOOpc scales of the parent simulation from which the 
initial conditions of this re-simulation were drawn, become signifi- 
cant only at the outer boundary of our re-simulation (see Appendix 
A2). Figure 7 shows that there are several sign changes in the torque 
profile, reflecting the specific state of the system at this time; the 



torque is time variable, but we find there are sign changes as well 
in a time-averaged sense. Overall, the net rate of change of the an- 
gular momentum of the gas very closely tracks the time-averaged 
gravitational torque from the stars. Hydrodynamic torques never 
induce very strong torques (greater than those shown), whereas the 
stellar gravitational torques can be a factor of ~ 10 — 100 larger 
than in Figure 7 at some times. 

The fact that the torques on the gas are dominated by stars 
is robust, and occurs for two reasons. First, the stars contribute 
significantly to the mass on the scales > 0. 1 pc that we focus on 
here (where star formation can occur). For typical star formation 
efficiencies, it is difficult to have the gas mass ^> 50% of the to- 
tal mass for a reasonable fraction of the lifetime of the system. On 
large scales, galaxies are known to not be so gas rich (although they 
certainly can reach ~ 50% gas fraction). On small scales, if the star 
formation efficiency is a few percent per dynamical time, then even 
a pure gas inflow from larger scales is likely to only remain gas- 
dominated for ~ 10 local dynamical times (on the smallest scales, 
~ 10 6 yr) - thus for the majority of the time during which inflow is 
continued, the system will contain a significant stellar mass. This is 
consistent with direct observations of sub-kpc regions of starburst 
galaxies (Downes & Solomon 1998; Bryant & Scoville 1999) and 
the ~ 1 — 10 pc nuclear scales around AGN (Hicks et al. 2009). 
On sufficiently small scales, < 0.01 — 0.1 pc, star formation will 
become inefficient, but at precisely those radii, we have (by defini- 
tion) essentially reached the a-disk. 

Even if the gas mass is comparable to or greater than the stellar 
mass, the self-torque of the gas on itself is much weaker than the 
torques from the stars on the gas. This has been demonstrated in 
detail in the case of large-scale perturbations from galaxy mergers, 
bars, and spiral waves (Noguchi 1988; Barnes & Hemquist 1996; 
Barnes 1998; Berentzen et al. 2007; Hopkins et al. 2009b). We 
show that it is also true for smaller radii in Figure 7. We discuss 
this in detail in a subsequent paper (Paper II; in preparation), but 
briefly outline the key physics here. 

It is well-known that the magnitude of the self-torque in a pure 
gaseous or pure stellar disk is second-order in the mode amplitude 
(oc \a\ 2 , where \a\ is the fractional magnitude of the asymmetry; 
see e.g. Lynden-Bell & Kalnajs 1972). Moreover, Kalnajs (1971) 
showed (with the exact solution for a logarithmic spiral) that the 
torques are further suppressed by powers of ~ kR and m, where 
k (m) is the radial (azimuthal) wavenumber of the mode. This is 
despite the fact that the instantaneous torque on gas/stars moving 
through a perturbation has no such suppression and is linear in \a\. 
For small mode amplitudes, the gas orbits periodically (in the test- 
particle limit) in response to the perturbation and the positive and 
negative linear terms exactly cancel because of the symmetry on 
either "side" of the mode. In a mixed stellar+gas system, however, 
the two perturbations are not exactly in phase because of dissipation 
in the gas. Perhaps most importantly, the gas can be driven into 
shocks by a sufficiently strong perturbation in the stars, leading to 
strong dissipation. As a result, the symmetric response of the gas 
is broken, and the net torque from stars on the gas is linear in \a\ 
and largest in the global limit. At the order of magnitude level, the 
inflow rate is thus given by 



M, 



\a\HgasR fi. 



(4) 



Note that expressed in terms of a viscosity v such that M ~ 
37r^E gas , equation (4) corresponds to v ~ (|a|/3-7r)V<J?, where V c = 
RSI is the local circular velocity. Equivalently, the inflow time of 
the gas is ~ ] « ] 1 1 . Typical values of \a\ are ~ 0.01 — 0.3 in our 
simulations when significant inflow is present (see §5). We stress 
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Figure 8. Properties of representative galactic (right panel), intermediate 
(middle), and nuclear-scale (left) simulations as a function of radius at a 
single time. The simulations are the same as in Figures 2 & 5 and the time 
chosen corresponds to near the peak of activity, with significant gas inflow 
to small radii. Note that the initial/boundary conditions for the smaller-scale 
simulations were taken from the larger-scale simulations shown, so the re- 
sults are reasonably continuous at the radial boundaries (despite the fact 
that these do not represent a self-consistent single simulation). Top: Gas 
and stellar surface density. Second-fmm-Top: Enclosed star formation rate 
in a given radius. Middle: Instantaneous inflow/outflow rate through the an- 
nulus (black line); note that there are several sign changes as a function of 
radius. Also shown is the order of magnitude estimate of the accretion rate 
produced by gravitational torques from a non-axisymmetric mode of am- 
plitude |a|, M ~ \a\ S gas i? 2 C (dashed red line; § 4). Second-from-Bottom: 
Scale height zo/R of the gas. Bottom: Toomre Q parameter of the gas. 



that, as described above, equation 4 differs fundamentally from the 
dimensional expectation for self-torques, which is second-order in 
\a\ and suppressed by terms in \kR\ and in. 

Figure 8 shows a number of properties of our fiducial sim- 
ulations from Figures 2 & 5, for each of the three spatial scales 
we consider: galactic (~ kpc), intermediate (~ 10 — 100 pc), and 
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Figure 9. Comparison of the true instantaneous inflow rate in the simula- 
tions (dM j&t ) to the simple dimensional scaling | a | Eg as R 2 £l(R) expected 
for pure gravitational torques, where \a\ is the magnitude of the strongest 
non-axisymmetric mode at a given radius R (directly measured in the simu- 
lations; with typical values ~ 0.01 — 0.3). Each solid line denotes a different 
simulation, with dotted lines showing the radii near the boundaries of our 
re-simulations, where the exact profile shape is suspect; colors correspond 
to initial gas fractions as in Figure 13. Thick dot-dashed lines correspond to 
ultra-high resolution simulations. The simulations are shown during the ac- 
tive/strong inflow phases, as in Figure 5. The inflow rate in the simulations 
is consistent with being produced by gravitational torques over 5 orders of 
magnitude in radius, with relatively little scatter (0.3 — 0.5 dex at all radii). 



nuclear (~ 0.1 — 10 pc). We show the stellar and gaseous surface 
densities, the local star formation rate, the gas inflow rate in com- 
parison to equation 4, and the vertical scale height and Toomre Q 
of the gas disk. We stress that Figure 8 shows three independent 
simulations, but that the smaller-scale simulations were initialized 
with properties drawn from the larger-scale simulations; this ex- 
plains the approximate, but not exact, continuity between the dif- 
ferent scales. Although Figure 8 is quite instructive, it can also be 
misleading to focus on the results of an individual simulation. The 
reason is that there is large variation in time and potentially large 
scatter introduced by modest differences in galaxy properties (due 
to, e.g., large-scale fragmentation of a spiral arm biasing the re- 
sults of a given simulation). For this reason we believe that it is 
extremely important to consider the statistical properties of a large 
number of simulations that vary the precise initial conditions, gas 
fraction, bulge to disk ratio, etc. In what follows we now discuss the 
statistical properties of the simulations summarized in Tables 1-3. 
In doing so, we will show a number of Figures that contain the re- 
sults of many of our simulations. In these Figures, the critical point 
to focus on is less the results of any given simulation (which can 
be hard to pick out), but rather the ensemble properties of all of 
the simulations (e.g., median, scatter, ...). Recall that the parameter 
space spanned by our simulations includes systematic surveys in 
the gas fraction, bulge (or BH)-to-disk ratio, and sub-grid equation 
of state of the gas. 

Figure 9 compares the order of magnitude estimate in equa- 
tion 4 with the true inflow rate in the simulations, as a function 
of radius; as in Figure 8, we show the results of simulations on 
all three of the spatial scales we have simulated, with dotted lines 
showing results near the radial boundary of a given calculation. In 
calculating M grav in equation 4, we approximate \a\ using the mag- 
nitude of the larger of the m = 1 and m = 2 Fourier component of 
the surface density distribution within R (see § 5). The results in 
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Figure 9 are at the peak of activity for each simulation with sig- 
nificant inflow, i.e., when M is near-maximum. Modulo a constant 
numerical prefactor, the simple scaling in equation (4) works well 
over the entire dynamic range of our simulations, with a scatter of 
~ 0.3 — 0.5 dex at any radius. The same equation also works for 
our ultra-high resolution simulations, which provide a smooth in- 
terpolation over the boundaries in our typical "re-simulation" cal- 
culations. The agreement in Figure 9 demonstrates clearly that the 
inflow rates do not follow the second-order expectations from self- 
torques in the absence of dissipation. In Paper II, we show that more 
detailed modeling of the gas response (accounting for e.g. the de- 
tails of the perturbation structure) can both predict the overall am- 
plitude of M in Figure 9 and reduce the scatter to < 0.3 dex. 

It is also important to note that the inflow rate in equation 4 
is independent of the sound speed c s , i.e., of the thermodynamics 
of the gas. Rather, the inflow rate is determined by the magni- 
tude of the non-axisymmetric gravitational potential. This is one 
of the reasons why we believe that our results are likely to be rel- 
atively robust, even given the significant uncertainties in the ISM 
physics. By contrast, if the transport were dominated by local vis- 
cous stresses, the accretion rate would be given by M v ; S e ~ 37r^E gas , 
where v ~ ac;/f2. For large-scale spiral waves, the angular mo- 
mentum flux associated with the wave corresponds to an effective 
viscosity that is at most v ~ c s r (Goodman 2003). Comparing this 
with the inflow rate induced by orbit crossing (eq. [4]) highlights 
that the dominant effect of the non-axisymmetric potential is that 
the stellar potential induces crossing orbits in the gas; the direct 
transport by spiral waves in the gas is smaller by at least a factor 
of ~ c s /V c . We find that in simulations that are initially 100% gas 
(even those initialized with already-strong mode amplitudes), the 
torques tend to remain fairly weak for a few dynamical times un- 
til the stellar component is non-negligible, at which point strong 
shocks appear, the stellar and gaseous modes develop a phase-shift, 
and the inflow rates rise sharply. The pure gas limit is, however, 
more sensitive to the thermodynamics of the gas and should be 
studied in more detail in future work. In our calculations, we find 
that it is critical to include stars, gas, and star formation to prop- 
erly capture inflow in galactic nuclei: it is the interplay between 
the stars and gas, and the competition between star formation and 
inflow, that sets the net accretion rate at small radii onto the BH. 

Our interpretation and analysis of the inflow rates assumes that 
the gas resides in a self-gravitating, at least modestly thin, disk. 
To demonstrate this explicitly, Figure 10 shows the vertical scale 
height and Toomre Q parameter of the gas as a function of radius. 
For each simulation, at a given radius R we determine the scale 
height zo of the gas by fitting the gas density distribution to a Gaus- 
sian with dispersion zo (after projecting to the plane perpendicular 
to the net angular momentum vector). At large radii, there is scatter 
owing to the difference between mergers (with sizeable zo /K) and 
unstable disks. From ~ 1 — 100 pc, zo /R ranges from ~ 0.05 — 0.25; 
this is determined by our assumed ISM equation of state, as well as 
kinematic features such as twists and misalignments. At the small- 
est radii za/R begins to decrease rapidly because the potential is 
increasingly dominated by the BH (the circular velocity increases 
but the maximum sound speeds are finite). Note that, at all radii, the 
minimum SPH smoothing length h sm \ < zo, so that the disk thick- 
ness is at least modestly resolved. Our simulated disks are always 
reasonably thin because the cooling time is always short compared 
to the dynamical time, so that the sound speed in the simulation 
reliably tracks that of our subgrid model. 

Figure 10 also shows the Toomre Q parameter of the gas, de- 
fined as Q = c s K/irGYid, where k is the local epicyclic frequency 
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Figure 10. Properties of the gas disks in our simulations. Top: Scale height 
Zo/R of the gas as a function of radius, for simulations with significant in- 
flows. Bottom: Toomre Q parameter of the gas. The inflows are robustly 
self-gravitating (Q ~ 1) and disky, with zq/R < 0.2. Each solid line denotes 
a different simulation, with dotted lines showing the radii near the bound- 
aries of our re-simulations, where the exact profile shape is suspect; colors 
correspond to initial gas fractions as in Figure 13. Thick lines correspond 
to ultra-high resolution simulations. The range in behavior emphasizes the 
importance of a large, statistically representative suite of simulations. 

and is the total surface density of the disk (gas plus stars). With 
some scatter, systems largely lie near Q ~ 1. This is broadly ex- 
pected for self-regulating disks; in the simulations, Q ~ 1 is a con- 
sequence of star formation, the assumed stellar feedback, and grav- 
itational dynamics. 

4.1 Gas Density Profiles 

Figure 1 1 shows the gas surface density as a function of time, at 
different radii, in the high-inflow simulation from Figures 2 & 5. 8 
Comparison of Figures 5 & 1 1 clearly indicates that there is a rea- 
sonable correlation between high inflow rates and high gas surface 
densities at the same radii. However, the correlation is not necessar- 
ily one-to-one because of the complex time-dependent dynamics on 
small scales. In some cases, the inflow leads the high surface den- 
sity (e.g., 0.4Myr at 0.1 pc or ~ 2 — 4Myr at lOpc), in which case 
the large M causes the high E, rather than the other way around. 
Occasionally, M and E can even be anti-correlated (e.g., at lOMyr 
on lOpc scales). 

Figure 12 shows the gas density profiles as a function of ra- 
dius for our ensemble of simulations (including the individual sim- 
ulations shown in Figures 5 & 6); for the gas-rich systems in the 
left panel, these surface density profiles are at the peak of activ- 
ity, i.e., when M is relatively high. During these active phases, the 

8 Note that for the re-simulations, the density at small r is intentionally 
initialized to be extremely small, but rapidly rises at very early times. 
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Figure 11. Gas surface density £ gas at different radii as a function of time, 
for the simulations shown in Figure 5. Top: Large scale galaxy merger 
simulation. The box shows the time interval re-simulated in the panel be- 
low. Middle: Intermediate-scale re-simulation. The box again shows the 
time interval for the smaller scale re- simulation. Bottom: Nuclear-scale re- 
simulation. Comparing with Figure 5, high inflow rates are generally corre- 
lated with higher local surface densities, but the relation is not one-to-one. 



gas reaches a quasi-equilibrium density distribution. For a given 
annulus, the net accretion rate to smaller radii is typically a small 
fraction of the total rate at which gas is supplied to that annulus. 
Thus to first approximation, the surface density can be estimated 
by considering the competition between star formation and inflow 
in a given annulus. If the star formation rate surface density scales 

3/2 

oc E gas , as in Kennicutt (1998), then setting the star formation rate 
inside some radius equal to the time-averaged inflow rate M m gives 
an expected (E gas ) oc M~{ 3 R~ 4 ^ 3 . Given M,„ oc R [ ^ 2 , which is rea- 
sonably consistent with Figure 5, this implies E gas oc R~ , similar 
to the results for the gas rich systems in Figures 11-12 (although 
slightly steeper than the numerical results). There is of course sig- 



nificant variation: we show the results from our large suite of simu- 
lations to emphasize the variation introduced by initial and bound- 
ary conditions, treatment of star formation and gas physics, etc., but 
also to highlight the robust average behavior and the dependence on 
global quantities such as gas fraction and inflow rate at larger radii 
(see Fig. 13 for the color scheme). 

Note that the gas surface densities achieved on small scales 
can be extremely large, ~ 10" — 10 12 Mq kpc -2 in the central 
~ 0.1 — lOpc; this is comparable to, or even somewhat larger than, 
the highest-surface density stellar systems known (Hopkins et al. 
2009e; Hopkins & Hernquist 2009). And it is a factor of ~ 10 4 
larger than the initial gas surface density at small radii (shown 
in Fig. 12 for comparison)! Assuming that a significant fraction 
of the gas eventually turns into stars, the relic stellar density pro- 
files expected from these gas-rich simulations are similar to the ob- 
served profiles of nearby "cusp" ellipticals (e.g., Kormendy et al. 
2009), which are indeed believed to be direct descendants of gas- 
rich mergers. However, estimating the "final" stellar profile at small 
radii, where our simulations are run for only a modest number of 
local dynamical times, requires a careful correction for the effects 
of duty cycle, so we defer a more detailed comparison to future 
work. 

Figure 12 (right panel) also shows the gas surface density pro- 
files for cases in which the initial large-scale gas inflows are not 
strong (e.g., Fig. 6). Most of these simulations represent either ex- 
tremely gas-poor mergers or (more commonly) systems that are 
weakly bar unstable or bar unstable but with significant bulge com- 
ponents (such that bars form efficiently, but inflow stalls at a large 
inner Linblad resonance). In such cases, the large-scale dynamics 
cannot increase the gas density at sub-kpc scales very far above the 
initial value of ~ 10 8 — 10 9 MQkpc~ 2 . Without this enhancement, 
the system is stable against secondary instabilities, and so there is 
little gas inflow at small radii (e.g., Fig. 6). 

In both the strong and weak inflow cases in Figure 12, our 
small subset of ultra-high resolution simulations (dot-dashed thick 
lines; 6 galaxy-scale runs that resolve down to ~ lOpc, and 5 
intermediate-scale runs that extend to < 1 pc) are fully consistent 
with the larger suite of lower resolution re-simulation runs. These 
high resolution simulations have continuous inflow from large to 
small radii, and can therefore be run self-consistently for a much 
larger number of dynamical times than our typical resimulation. 
In spite of this advantage, we find that our resimulation technique 
yields very similar results. 

5 CONDITIONS FOR (IN)STABILITY 

In the preceding sections, we have presented results for cases in 
which gas inflow from large scales both is, and is not, sufficient 
to trigger secondary instabilities at small radii, leading to efficient 
gas inflow down to ~ 0. 1 pc. Here we provide a more quantitative 
assessment of the conditions under which there is significant inflow 
of gas to sub-pc scales. 

Figure 13 shows how several different measures of the effi- 
ciency of angular momentum transport and gas inflow depend on 
the initial gas fraction / gas and the "bulge" to total ratio B/T - the 
same parameters that strongly influence the morphology of the gas 
(Fig. 3 & 4). We show results for a few different equations of state 
^eos (see Fig. 1). The quantities we use to define the efficiency of gas 
inflow are the gas mass within 10 and 0.1 pc M gas (< R) (top panel), 
the inflow rates within the same radii (middle panel) and the for- 
mal amplitude of the non-axisymmetric gravitational perturbations 
at ~ 1 and ~ 100 pc (bottom panel). For the inflow rate (middle 
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R [pc] R [pc] 

Figure 12. Left: Gas surface density profiles during the strong inflow phases, i.e., when the accretion rate at small radii peaks. Each solid line denotes a 
different simulation, with dotted lines showing the radii near the boundaries of our re-simulations, where the exact profile shape is suspect; colors correspond 
to initial gas fractions as in Figure 13. Thick lines are the ultra-high resolution simulations. Each simulation consists of re-simulations of the nuclear dynamics 
during gas-rich galaxy-galaxy mergers, as in Figures 5 & 11. Dashed red line shows the initial gas density profile of the large-scale simulations (run b3ex(ic) 
in Table 1). The gas density profile is quasi-steady over the entire active phase, with star formation offsetting inflow; the time variation in S gas for one example 
is shown in Figure 11. Right: Gas surface density profiles for simulations of isolated barred-disk galaxies and the corresponding re-simulations. In cases with 
a very small pre-existing bulge, the central gas density can increase by an order of magnitude, but not the 2 — 4 orders of magnitude seen in the left panel. 



panel), we show both the time average and "peak" values to convey 
the variability as a function of time in the simulation. 

The mode amplitude is measured in the standard fashion from 
the surface density distribution at a given radius, using 



\a m (R,t) \ = 



||/ o 2 "S(/?,0) exp(/m^)d0|| 



(5) 



We measure the amplitude from the stellar disk since it is slightly 
more robust to local clumping, but using the gas surface density 
gives similar results. For simplicity, we show the results for the 
largest-amplitude mode in each case, but this is almost always m = 
1 in the nuclear-scale simulations (right panels) and an m — 2 bar 
or m — 1 spiral in the intermediate-scale simulations (left panels). 
We measure the relevant amplitude at radii slightly larger than the 
radii where the inflow is measured, since it is this non-axisymmetry 
that drives the inflow (it is also somewhat less noisy at larger radii). 

The "bulge-to-total" ratio in Figure 13 is defined as the mass 
fraction in all spherical components within the specified radius, 
since these all serve to suppress instability in the same manner (e.g. 
any black hole, bulge, halo, and/or nuclear star cluster would qual- 
ify). Gas fraction is defined as the gas fraction of the disky compo- 
nent within the same radius. 

Figure 13 shows clearly that there is a very strong trend of 
inflow strength with B/T, in the expected sense. Strongly disk- 
dominated systems have characteristic inflow rates 3 orders of mag- 
nitude higher than strongly "bulge-dominated" systems; likewise 
mode amplitudes and gas masses inside small radii are higher by 
1—2 orders of magnitude. 

There are several interesting additional results in Figure 13. 
First, within the range of g cos that we consider (motivated by the 
observations in Figure 1), the inflow properties do not depend sig- 
nificantly on the ISM physics. This is consistent with our previous 
arguments; so long as the sub-resolution model is sufficient to pre- 
vent catastrophic fragmentation and runaway star formation, gravi- 



tational torques that are reasonably independent of the ISM micro- 
physics dominate the angular momentum redistribution. Note that 
there is also no distinguishable difference between the results of 
our ultra-high resolution runs (shown with large symbols) and our 
standard resimulations. Second, the inflow properties depend much 
more strongly on B/T than on the initial gas fraction; this is because 
the stellar disk generates the same instabilities as the gas disk. In 
fact, the collisional nature of gas means it cannot support certain 
strong modes (Toomre 1977), so the transport can be less efficient 
in the fully gas-dominated limit (the gas-rich systems modeled here 
often do not develop strong modes until a non-trivial fraction of the 
gas has turned into stars). 

Third, there is considerable scatter at all B/T. This is in no 
small part due to the time-variable nature of the small scale dynam- 
ics (§3). In addition, however, a quantity such as B/T is not a global 
invariant; it depends on both position and time in these simulations. 
As such, some simulations with a given B/T at the specified radius 
can have a larger or smaller B/T at some larger radius, leading to 
more or less efficient inflow from those scales that in turn leads to 
more or less efficient inflow to smaller scales. Identifying B/T at 
a given radius and time is thus a crude measure of stability. More- 
over, the torques and inflow rates depend on the precise structure of 
the modes involved; a significant portion of the scatter in Figure 13 
reflects the difference between e.g. bars, spiral arms, and loosely or 
tightly-wrapped waves. This large scatter highlights the importance 
of using a large suite of simulations to conduct parameter surveys 
and build statistical samples - any individual system, studied for a 
small amount of time, can be highly non-representative. 

5.1 Intermediate (sub-kpc) Scales 

For the intermediate scale results shown in the left panel of Fig- 
ure 13, we can gain some analytic understanding of the behavior 
by following Shlosman et al. (1989), who discuss the criteria for 
secondary instabilities ("bars within bars"). It is well established 
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Figure 13. Inflow strengths as a function of bulge to total ratio (B/T) for several initial gas fractions (colors) and ISM equations of state g C os (symbol type), 
for a subset of our simulations. Large symbols in each refer to the subset of ultra-high resolution simulations. Top: Total gas mass inside an inner radius of 
lOpc (left; from our intermediate-scale re-simulations) or 0. 1 pc (right; from our nuclear-scale re-simulations), near the peak of accretion in each simulation. 
"Bulge" (B) refers to any spherical term in the potential (bulge+halo+black hole). Middle: Gas inflow rate through the relevant inner radius. Solid points denote 
the time-averaged inflow rate over the full re-simulation while open points are the maximum seen in any individual snapshot (to give some sense of the duty 
cycle). Bottom: Amplitude of the dominant non-axisymmetric mode in the stars (generally m = 2 at intermediate scales, m = 1 at nuclear scales), averaged 
over the peak/strong inflow phase. Inflow properties scale broadly with B/T as expected, but the large scatter and sensitivity to other parameters makes a large 
parameter survey critical. 



that a self-gravitating disk becomes unstable to large-scale gravi- 
tational modes roughly when the Ostriker-Peebles criterion is met, 
i.e., when 7rot/|W| > tan, where T IBt ~ M d V} is the kinetic energy 
of rotation (in the rotationally supported - i.e. kinematically cold 
- component), \W\ ~ GM 2 /R is the potential energy, and lent is a 
threshold value ~ 0.15 — 0.26, depending on the nature of the in- 
stability (e.g., Bardeen 1975; Aoki et al. 1979). If a mass fraction 
fd — 1 — B/T is in the disk, the condition for instability generically 



becomes 

fd>fan(taK,a/R) (6) 

where a is the bulge scale length, R is the disk scale length, and the 
exact functional form of / c ,i t depends on the profiles of the disk and 
bulge. For Kuzmin and Plummer profiles for the disk and bulge, 
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respectively, we find that this instability criterion becomes 



M d [<R] 



> 



M h +M d [<R] \4[l-2t cnt ]a 



3nt c , 



R 



1/2 



16 a 



(*«a) (7) 
(if- a) (8) 



where M d [< R] is the total (stars + gas) disk mass within R and M h 
is the total bulge mass. 

If the disky material of interest has a scale-length compara- 
ble to that of the bulge (generally true for the intermediate- scale 
simulations here), equation (6) simply becomes B/T < 1 — r cr ii ~ 

0. 7 — 0.8. This is reasonably consistent with the fact that the inflow 
rate and interior gas mass in the intermediate scale simulations in 
Figure 13 (left column) asymptote for B/T > 0.6 — 0.7. 

Shlosman et al. (1989) show that an initial gas disk must be 
compressed to at least R wvl /Ri ~ (1.0 - 2.5)/ g 2 as /[l + 0.2 (B/T) 2 ] 
to satisfy equation (6), again for systems in which the disk and 
bulge initially have similar scale-lengths. For the case of major 
mergers, we can compare this to Covington et al. (2008) and Hop- 
kins et al. (2009a,c)'s estimates of how much angular momentum 
the gas loses during the merger. These authors find that the gas of- 
ten inflows until it becomes self-gravitating - as a result, secondary 
bar instabilities should be ubiquitous. More quantitatively, they find 
Rww/Ri « /gas/(l +/gas//o) with / « 0.2-0.3. Using this, we es- 
timate that mergers with reasonable gas fractions > 0.3 — 0.5 will 
lead to secondary instabilities and rapid inflow. For isolated disk 
galaxies with bars, the critical criterion is whether the inner Lin- 
blad resonance lies inside or outside of R new ; this will generally not 
be true in low-/ gas or high-B/T systems. 

5.2 Nuclear (pc) Scales 

Inside the radius where the BH begins to dominate the gravita- 
tional potential, it becomes increasingly difficult for system to be 
globally gravitationally unstable in the sense of r rot /|VF| > f crit . 
Specifically, in the potential of the BH this criterion implies that 
the disk would no longer be unstable inside a radius R m i„ ~ 
lOpc(M BH /lO 8 M ) 1/2 (E rf /lO 11 M kpc- 2 )- 1/2 . This confirms 
our earlier claims (e.g., §3.3) that the character of the instabilities 
responsible for angular momentum transport must change near the 
BHs radius of influence. Indeed, we see numerically that the trans- 
port on small scales is dominated by m = 1 modes (e.g., Fig. 4). 

A full discussion of the origin of the m — 1 modes is beyond 
the scope of this paper; we will present a detailed analytic analy- 
sis (in simplified BH+disk models) of their structure, growth rates, 
and pattern speeds in future work (in preparation). We can, how- 
ever, present a simple analysis indicating how we believe the modes 
arise. For m = 1 modes, — n/m— > in a Keplerian potential. It is 
well-known that this allows very low frequency (low pattern speed 
fl p ) modes to be present in a Keplerian potential (e.g., Tremaine 
2001). Moreover, Tremaine (2001) showed that the only modes that 
can be global and exert coherent strong torques over a large dy- 
namic range in a quasi-Keplerian potential are these m = 1 modes. 
We find that the modes in the simulations indeed have small Q,,, 

1. e., Q p <§C Q and m = 1, deep in the potential well of the BH. Note 
that for a collisionless particle, rather than the stellar+fluid disks of 
interest here, this simply corresponds to a standing elliptical Ke- 
plerian orbit. Small fi^/fi (and wavenumber m — 1) is important 
because it implies near-resonance conditions at small radii; this al- 
lows the stars to produce strong torques on the gas, which drives or- 
bit crossings, shocks, and angular momentum loss, even when the 
mass of the disk at small radii is much less than the BH mass (see 
e.g. Chang et al. 2007). Following Tremaine (2001) and expanding 



the equations of motion for a linear mode in a nearly-Keplerian po- 
tential (in the WKB limit) it is straightforward to show that the 
magnitude of the induced velocities (8v/V c ) from a given mode 
scale as 5v/V c ~ (<5E/E)M rf (< R)/M m = \a\M d (< R)/M B h for 
all m =fc 1; form = 1, however, 8v/V c ~ \a\. In other words, the con- 
ventional wisdom that torques from global modes are strongly sup- 
pressed in the spherical potential of a BH is generally true, but the 
resonance condition that k ~ Q. in Keplerian potentials allows for 
large coherent eccentricities, self-gravitating collective effects, and 
torques from m = 1 modes at arbitrarily small disk-to-BH masses. 
As a result, although we certainly see higher-m modes in our small- 
scale simulations, the global structure is almost always dominated 
by a coherent global m = 1 eccentric disk. 

However, Tremaine (2001) also showed that slow modes with 
pattern speeds Q. p <§C Q,, are linearly stable in quasi-Keplerian po- 
tentials. How, then, do these modes arise in the simulations? To 
start, we consider the dispersion relation for linear waves in the 
WKB limit (see e.g. Lau & Bertin 1978), which is given by 



(^) 2 = l + ,-2/M(|^W)%A(W + 



2/ 2 
777 7 7 



where 



ttE,;/? 2 _ M d (<R) 



M mc (<R) M tot (<i?) : 



A=l + 



2ra 2 (3- 



{f+l)Qkr\ 2 + m 2 )' 



1+2 



d\nV c dhiM enc (</f) 



(9) 
(10) 

(11) 
(12) 



d\nR dlnfl 

and h — c s /Q is the disk thickness. We now consider the limit of 
equation (9) such that kr is modest. Of course, our analytic treat- 
ment of such quasi-global gravitational modes should be consid- 
ered with due caution, as the WKB approximation necessary to de- 
rive the dispersion relation is not justified. But this nevertheless 
points towards the conditions under which global gravitational in- 
stability may be present. 9 In this limit, equation (9) will have grow- 
ing modes if the right-hand side is negative, which for a cold disk 
(h <C r) requires 



h> 



1 (1+7,) 2 



(13) 



2m 7 — 7/ 

As v — > 0, i.e., as the potential becomes Keplerian, equation 13 
suggests that m = 1 modes can become unstable when f d > 0.1. 
This is consistent with the fact that we find the m — 1 modes at 
nearly all B/T (Fig. 13). Equation 13 also allows for the instability 
of higher-m modes, but these cannot exert coherent strong torques 
on the gas, nor can they propagate efficiently to small radii, so they 
do not dominate the structure or inflows we see in our simulations. 

Because the disk in general will have more mass at large radii, 
any instability will grow outside-in; for massive disks, the most- 
unstable point is where M d (< R) ~ 0.5 - 1M BH (fd ~ 1/3 - 1/2). 

9 Note also that equation 9 is derived with the WKB approximation, but in- 
cluding terms up to second order in |tcR| _1 (dropping out-of-phase terms), 
which makes it somewhat more accurate for the large-scale modes of inter- 
est here. Specifically, this leads to the additional enhancement of the growth 
rate A ~ 1 + T sin;', where T = 8\nQ/d\nR and i is the arm pitch angle; 
this approximates (at this order) the effects of the swing amplifier, an im- 
portant contribution to the instability criterion. 
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Here modes can grow on a dynamical time. Equation 9 shows that 
this growth is for fast modes: overstabilities with Re(o;) ~ mfl. In 
the simulations, the slow modes typically dominate the torques at 
sub-pc scales. However, we find that the fast mode itself drives 
the slow modes. Specifically, the m = 1 mode first appears around 
the most unstable point ^ cr i t predicted above, where Mj(< R cra ) ~ 
0.5 — IMbh- This rapidly leads to an eccentric disk at these radii. 
As the instability causes gas to lose angular momentum, the gas 
falls inwards, and begins turning into stars. The eccentric mass dis- 
tribution at i? C rit then in turn efficiently induces a similar eccen- 
tricity at smaller radii. As noted in Tremaine (2001), the m — 1 
slow mode may be linearly stable, but it can be excited in the first 
place - and will be long-lived - in response to an external asym- 
metry in the potential. Our interpretation of the simulations is that 
the unstable self-gravitating disk near ~ R C n t provides the asym- 
metry that generates the slow, low pattern speed, m = 1 mode at 
smaller radii. The pattern speed £l p — f2(/f C rit) is conserved as the 
mode moves in, but the circular speed Q(R) is rapidly increasing 
at smaller radii near the BH, so the mode at small radii is indeed 
a slow mode - for the typical parameters here, the pattern speed is 
fi p ~ 1 — 5 km s~ 1 pc~ 1 , set by the circular speeds at ~ 10 — lOOpc, 
where Mj(< R) ~ Mbh (for comparison, Cl > 1000km s~ 1 pc~' at 
R < 1 pc). The slow mode at small radii responsible for the angular 
momentum transport is therefore not a local small-scale instabil- 
ity, but part of a global instability beginning at larger radii. In the 
future, we plan to explore the origin of this global m — 1 mode 
in more detail. Calculating the dynamics of the gas and stars, and 
including self-gravity and star formation, are all important for cap- 
turing the inward propagation of these modes. 

Similar m — 1 modes have been studied previously as a mech- 
anism for accretion in fluid disks (Adams et al. 1989; Shu et al. 
1990; Ostriker et al. 1992; Bournaud et al. 2005), both in the galac- 
tic context and in proto-planetary or proto-stellar disks (e.g., Baner- 
jee et al. 2004; Boley et al. 2006; Vorobyov & Basu 2006, 2009; 
Krumholz et al. 2007; Cai et al. 2008). In a purely gaseous disk, 
wave packets can propagate through the OLR to r — > oo, eventually 
becoming simple sound waves (Adams et al. 1989). Because the 
waves can freely escape to the outer boundary, infinite pure gaseous 
disks in nearly-Keplerian potentials do not support strong growing 
m — 1 modes, and the mode growth is quite sensitive to the descrip- 
tion of the disk boundary. However, for a stellar disk, the mode can- 
not propagate beyond the OLR where A = k 2 — m (VL — Q.,,) 2 = 0. 
Refraction of the stellar waves at this boundary is important, and 
implies that global m — 1 modes can grow even when the disk ex- 
tends to R 3> Rolr- As a result, we shall show in subsequent work 
that the m = 1 modes in quasi-Keplerian stellar disks are relatively 
insensitive to the outer boundary conditions. This is also implicit 
in Tremaine (2001), where for many of the modes the outer radius 
of the disk could be arbitrarily large. The ability of stellar m = 1 
modes to grow more readily again highlights the importance of con- 
sidering the two component stellar + gas system when considering 
gas inflow in galactic nuclei. 



6 ACCRETION VERSUS STAR FORMATION 

Given the gas inflows modeled here, our results imply a correlation 
between the accretion rate onto the central BH - and hence the 
AGN luminosity - and the star formation rates on larger scales. 
Figure 14 shows this correlation for a number of our simulations. 
We take the inflow rate at 0.1 pc as the AGN accretion rate, which 



corresponds to a bolometric luminosity of 



Z*oi,agn = e r M B HC « 5.7 x 10 



M BH 
Mq yr" 1 



ergs 



(14) 

where the radiative efficiency is e r = 0.1. Figure 14 compares the 
BH accretion rate to the total star formation rate inside a (pro- 
jected) radius R, at a few representative radii: true nuclear scales 
(R < lOpc), more readily observable nuclear scales (R < lOOpc), 
the smallest scales observable in most moderate to high-redshift 
systems (R < lkpc), and the total star formation rate, integrated 
over all radii (R — > oo). For each point, in Figure 14, the projected 
SFR(< R) is technically the median over ~ 100 sightlines, but pro- 
jection effects are negligible. For the larger radii, we measure the 
SFR in an appropriate intermediate or galaxy-scale simulation. The 
BH accretion rate is then determined from that in the appropriate 
re-simulation or set of re-simulations matched to the galaxy scale 
simulation central properties at the given time. Because a given dy- 
namical time in a galaxy-scale simulation corresponds to many dy- 
namical times (and simulation outputs) on very small scales, the 
number of points is larger in the plots going out to larger scale. 

Not surprisingly, Figure 14 shows that there is a correlation 
between accretion rate and star formation rate, on all scales. There 
is significant scatter in this correlation, although the scatter is less 
when the SFR is evaluated inside a small nuclear annulus. If we 
assume a linear proportionality between BH accretion rate and star 
formation rate in some annulus, we find 



M B H~M«(tf < lOpc) 

~0.1M*(fl < lOOpc) 
~0.03M»(tf < lkpc) 
~ 0.003 M, (total) 



(15) 
(16) 
(17) 
(18) 



for the median relation in the simulations, which corre- 
sponds to a bolometric luminosity ratio of Lboi,SF/£boi,AGN ~ 
(0.01,0.1,0.3,2.5) at (R < 10pc,/{ < 100pc,^ < lkpc,tf < 
oo), respectively, or a ratio of the star formation-powered in- 
frared luminosity to the hard X-ray luminosity from the AGN of 
£ir,sf/£hx,agn ~ (0.2, 2, 6, 60) in the same annuli. 

In greater detail, the correlations are not precisely linear. Pa- 
rameterizing the correlation between star formation and BH ac- 
cretion as a power-law, Mbh oc M™, we find that the correlations 
above are approximately valid for Mbh ~ 0. 1 — IOMq yr~ 1 , but a 
changes from near-linear at small scales (R < lOpc) to somewhat 
super-linear (a ~ 1.5) on large scales. For the total star formation 
rate, we find 



MBH~O.5M yr->( 5o]i ^ vFT ) 



5OM yr- 



(19) 



However, we caution that our simulations do not include the ap- 
propriate physics to resolve the variety of processes that can pro- 
duce low accretion rates <§C 0. 1 Mq yr~ 1 , nor do they sample a large 
regime of galaxies with total star formation rates < 1 Mq yr~ 1 , so 
the fits above should not be used in the limits of weak accretion 
and star formation. The qualitative point is nonetheless important: 
non-axisymmetric instabilities due to self-gravity become ineffi- 
cient at low gas densities (low star formation rates), and so other 
mechanisms must be important for fueling at least some of the low- 
luminosity AGN population. 

It is notable that even during the active phases modeled here, 
the average relation between accretion rate and total SFR is such 
that the systems are most often not AGN dominated in a bolometric 
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Figure 14. The predicted relationship between the BH accretion rate (BHAR) and the star formation rate (SFR) inside a given radius. For each simulation, 
we plot the SFR integrated within the indicated projected radius, at different times during the simulation, compared with the < 0.1 pc BH accretion rate 
from appropriate re-simulation. The characteristic timescale for variability at small radii - which sets the variability in the BHAR - is much shorter than the 
corresponding timescale at large radii; at fixed SFR, each system can thus have a variety of BHARs. Nevertheless, the BHAR correlates with SFR, with the 
least scatter for the nuclear star formation (upper left). The correlations seen here are a natural consequence of the fact that high gas densities, and thus high 
star formation rates, are necessary to trigger secondary instabilities and drive accretion on small scales. 



sense (although the two are certainly comparable, in these phases). 
In other words, even the most active systems will be consistent with 
SFR-dominated properties except for a small fraction of the time; 
this is consistent with observations such as the FIR-radio correla- 
tion (Walter et al. 2009; Riechers et al. 2009; Wang et al. 2008). 
However, the scatter in accretion rates at fixed SFR is large, and 
most of the growth in BH mass comes from the rare, short-lived 
times when the systems scatter to higher Lagn/£-sf- Roughly, we 
can estimate these duty cycles with the scatter in Figure 14 - about 
~ 30% of the systems in "peak" AGN phases have sufficient accre- 
tion rates (M BH > 0.01 M,) so as to be AGN-dominated. 

The predictions in Figure 14 can be compared to a number of 
observations. Wu & Cao (2006) compile total star formation rate 
estimators in ~ 100 nearby AGN ranging in luminosity from LIN- 
ERS (Mbh < 0.01 M yr~ 1 ) to PG quasars (M B h ~ 1 - 5M Q yr~ 1 ). 
The sample details are discussed in Satyapal et al. (2005), but the 
authors find they are reasonably robust to the choice of SFR or BH 
accretion rate indicator. Silverman et al. (2008) consider a simi- 
lar compilation in more luminous systems in the COSMOS sur- 
vey, spanning a redshift range z ~ 0.3 — 1.1. At total SFR values 
> lOMQyr" 1 , these samples agree fairly well. Specifically, for 
M„ ~ 10 — 1000 Mq yr~ 1 , these different compilations are consis- 
tent with a median Mbh ~ 0.001 — 0.01 M*, with factor of ~ 3 dis- 
persion; this is quite similar to our results shown in Figure 14. On 



the other hand, Ho (2005) used Oil line emission to argue that the 
star formation rate associated with luminous type 1 AGN is at best 
comparable to the BH accretion rate, inconsistent with the median 
integrated galaxy predictions in Figure 14. This may be because ob- 
scuration makes Oil an imperfect proxy for star formation in dense, 
obscured, starbursts. In addition, luminous Type 1 AGN are prob- 
ably, by selection, on the tail end of the BHAR/SFR distribution 
(and could in principle be in an AGN feedback-dominated stage of 
evolution not modeled here). Indeed, Kim et al. (2006) followed up 
Ho's work and found enhanced star formation in type 2 quasars, 
similar to our predictions here. 

Wang et al. (2007) consider spatially resolved estimates of 
nuclear luminosity and emission-line (SFR) profiles in a range of 
Seyferts and quasars (mostly at low redshift) at spatial scales from 
~ 0.05 - 3 kpc (see Imanishi 2002, 2003; Imanishi & Wada 2004). 
Their results are similar to those in a variety of other studies (see 
e.g. Kauffmann et al. 2003; Kewley et al. 2006; Evans et al. 2006; 
Schawinski et al. 2009). At both R < lOOpc (where the observa- 
tions span a dynamic range in SFR(< 100 pc) ~ 0.1 — 10Mq yr _1 ) 
and/?< lkpc (SFR(< lkpc)~ 1 - lOOM yr~'), our results are 
reasonably consistent with these observations, in both normaliza- 
tion and scatter. A small sample is studied at extremely high res- 
olution in Davies et al. (2007), where star formation and AGN lu- 
minosity are resolved at ~ 10 — 30pc scales. They find that, within 
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these radii, the luminosity from star formation is ~ 1 — 5% of that 
from the BH, again consistent with our prediction here. 

7 DISCUSSION 

We have presented high resolution smoothed particle hydrody- 
namic simulations of the dynamics of gas in the central regions of 
galaxies, in order to study gas inflow from galactic scales ~ lOkpc 
down to the scales of a canonical AGN accretion disk (~ O.lpc). 
We use the properties of galactic-scale simulations - both merg- 
ers and isolated galaxy disks - to initialize new higher resolution 
"re-simulations" that focus on the < 100 pc dynamics. This tech- 
nique is then applied a second time to follow the gas to ~ 0. 1 pc. 
By intention, our calculations do not employ particle splitting and 
are not all exact realizations of the small-scale dynamics associated 
with a single, specific larger-scale simulation. Our "re-simulations" 
should not be taken as such an exact realization. Instead, since it is 
not practical to directly solve the "full" problem (from ~ lOkpc to 
10 _5 pc), we have focused on understanding and isolating the key 
physics involved. To provide a more physical model, however, we 
use initial conditions drawn from a variety of larger galaxy-scale 
simulations, rather than an ad hoc set-up. This approach allows 
us to systematically survey, for the first time, how a large variety 
of initial conditions and galaxy properties affect gas dynamics in 
galactic nuclei. 

It is still not feasible, with any numerical scheme, to simul- 
taneously model the small scales of an AGN accretion disk and 
those of a galaxy for cosmological timescales. For this reason, di- 
rect "zoom-in" approaches have thus far been fairly limited, usually 
evolving each region for just a few local dynamical times. Our nu- 
merical approach attempts to overcome this severe limitation. First, 
the smallest regions of each simulation are simulated for many lo- 
cal dynamical times. Even more importantly, we quantity the dy- 
namics on small scales using ~ 100 simulations of a large range 
of initial conditions and galaxy properties This approach does not 
require an exact mapping between one large-scale simulation and 
another on smaller scales (although we often interpret our results in 
this manner). For example, even if our nuclear-scale simulations of 
the central ~ lOOpc of a galaxy are not self-consistently matched to 
those of a specific larger-scale simulation, they still represent valid 
simulations of a given set of initial conditions of a bulge+stellar 
disk+gas disk+BH system, evolved for many dynamical times. As 
such, we can still use them to understand the physics of angular 
momentum transport and BH fueling at these radii (given these ini- 
tial conditions). This methodology allows us to survey a wide range 
of initial conditions and galaxy properties and isolate the physical 
parameters that are most important for the physics of gas inflows 
on small scales in galactic nuclei. We have also checked our "re- 
simulation" approach by carrying out a small number of ultra-high 
resolution simulations that bridge the gap between different resim- 
ulated regions; we find in all cases that these ultra-high resolution 
simulations validate our resimulation methodology. 

Our calculations demonstrate that for sufficiently gas-rich, 
disk-dominated systems, which have a large inflow of gas into 
~kpc scales, there is a cascade of non-axisymmetric gravitational 
instabilities that ultimately produces BH accretion rates as high as 
~ lOMQyr -1 at < 0.1 pc (Fig. 13). Moreover, we have explicitly 
shown that these conditions are satisfied during major mergers of 
gas-rich galaxies, thus providing theoretical support for the con- 
nection between mergers and BH growth. It is, however, also im- 
portant to stress that our work only implies that galaxy mergers are 
a sufficient condition for fueling a central quasar, not that they are 



necessary: similar inflows can be obtained in isolated gas-rich, bar 
or spiral wave-unstable galactic disks. 

In broad terms, our results support the "bars within bars" sce- 
nario proposed in Shlosman et al. (1989). However, we show that 
the secondary instabilities on intermediate scales (~ 10— 100 pc) 
exhibit a diverse range of morphologies, including standard nuclear 
spirals, bars, rings, barred rings, one or three-armed spirals, and 
irregular/clumpy disks and streams (Fig. 2 & 3). This is very im- 
portant for comparing to real galaxies: observations have generally 
found that nuclear bars are not ubiquitous in AGN (perhaps not 
substantially more common than in non-active systems). However, 
these same observations have found that some form of asymmet- 
ric nuclear structure is ubiquitous, with the most common features 
being nuclear spirals and rings very similar to those seen in Fig- 
ure 3 (see Martini & Pogge 1999; Martini et al. 2003; Peletier et al. 
1999; Knapen et al. 2000; Laine et al. 2002; Knapen et al. 2002; 
Greene et al. 2009, and references therein). Our work shows that 
these observations are not necessarily in conflict with the hypothe- 
sis that gravitational torques dominate inflow on sub-kpc scales in 
AGN. Rather there is a large diversity in the observational mani- 
festation of these gravitational torques. "Bars within bars" should 
not be taken too literally. Instead a more accurate characterization 
would be: "it's non-axisymmetric features all the way down" (or 
"stuff within stuff"). 

This diversity in observational appearance owes both to the 
effects of global parameters such as gas fraction and bulge-to- 
total ratio, and more subtle variations induced by the relative scale 
lengths of disks and bulges, the precise profile shapes, and the 
thickness/dispersion in the sub-kpc bulge, stellar disk, and gaseous 
disk components. Observationally, these parameters are all seen to 
vary considerably even in isolated galaxies, let alone in chaotic 
merging systems; as a result, we expect the morphological diver- 
sity exhibited in Figures 2 & 3 to be the norm. 

Once the gas reaches the BH radius of influence (~ lOpc for 
a typical ~ L, system), the potential becomes quasi-Keplerian and 
spherical (dominated by the BH itself); the system is thus no longer 
bar unstable. However, the gas is still locally self-gravitating and 
prone to forming stars if it does not inflow rapidly. Indeed, this is 
traditionally the range of radii at which it has been the most difficult 
to produce the large accretion rates needed for luminous AGN and 
quasars (Shlosman et al. 1990; Thompson et al. 2005). We find, 
however, that new large-scale non-axisymmetric modes arise at ~ 
pc-scales and robustly allow for continued gas transport to smaller 
radii. Specifically, for sufficiently gas-rich, disky systems, we find 
an eccentric/lopsided disk or a one-armed spiral instability (an m = 
1 mode); see Figure 4. This leads to large inflow rates of ~ 1 — 
lOMQyr -1 into the central 0.1 pc, comparable to what is needed 
to explain the brightest quasars observed. At radii < 0.01 — 0. 1 pc, 
the accretion flow is no longer self-gravitating and local viscous 
heating is sufficient to maintain Q > 1 and transport gas to the BH 
(Goodman 2003; Thompson et al. 2005); the disk is approaching a 
canonical thin Keplerian accretion disk. 

The m — 1 modes seen here have been studied previously as 
a mechanism for accretion onto BHs and protostars (Adams et al. 
1989; Shu et al. 1990; Ostriker et al. 1992), but have largely been 
neglected in studies of fueling luminous AGN. This is in part be- 
cause most previous simulations neglected either star formation or 
the self-gravity of the gas and stars (owing to computational limits 
or, of course, the fact that these complications are not present in the 
case of a protostar). Moreover, previous calculations have shown 
that although near-Keplerian potentials can support the "slow" pat- 
tern speed m — 1 modes we find here, such modes are linearly 
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stable (Tremaine 2001). Why, then, are these modes ubiquitous in 
our simulations? We find that the inclusion of stars and gas, self- 
gravity, live star formation, and some model for stellar feedback are 
all important for the behavior of these instabilities. In particular, we 
typically find a complex interplay between the gas and stars (§4 & 
5.2). The system often develops the m — 1 modes first near co- 
rotation, where the disk and BH contribute comparably to the po- 
tential; the mode growth is faster in the stellar distribution, which 
can support self-crossing orbits and which is much less sensitive 
to the outer properties of the disk. As gas streams move through 
the lopsided stellar distribution, they are torqued, experience or- 
bit crossing, and lose angular momentum and energy (as in, e.g., 
Chang et al. 2007). As the gas moves to smaller radii, the stars at 
larger radii provide a less efficient angular momentum sink. How- 
ever, star formation ensures that new stars are formed in situ (them- 
selves in a lopsided distribution), allowing the gas to continue to 
flow in. This process leads to the m = 1 mode propagating inwards 
into the gravitational potential well of the BH, from the larger radii 
where it originated. The torques and inflow rates are thus quite dif- 
ferent in stellar+gas systems as compared to pure gas disks. In par- 
ticular, orbit crossing induced by the stars can lead to much more 
efficient gas inflow; we will study this analytically in Paper II. 

Our calculations show that the inflow of gas on sub-kpc scales, 
including the BH accretion rate at < 0. 1 pc, is highly time vari- 
able (Fig. 5). This variability is a result of gravitationally unstable 
modes forming, damping, and, re-forming, and the fact that individ- 
ual clumps or dense regions can dominate the inflow rate at a given 
time. The "duty cycle" for large inflow rates is modest, even in sys- 
tems with a large time averaged accretion rate. In one of our most 
violent simulations, a massive nuclear ring forms at ~kpc; owing 
to an inner Linblad-like resonance, the system has an outflow from 
sub-kpc scales, and inflow from the merger at larger radii. Since 
the outflow/inflow rates are large, the ring soon becomes strongly 
self-gravitating and globally collapses, leading to one of the largest 
net accretion rates into ~ lOpc in any of our simulations. How- 
ever, more than ~ 95% of the time during the "active" phase, the 
system has strong outflow from the sub-kpc region. Moreover, be- 
cause the gravitational modes driving accretion at large scales (e.g. 
the merger), at ~ 100 pc (secondary instabilities), and at ~pc (the 
lopsided nuclear disk) are physically separate, their variability on 
short timescales is not fully coupled. This is clearly important when 
relating observed torques on large scales to AGN accretion rates on 
small scales. In many observed AGN or barred disks (with a proper 
inner Linblad resonance), there are no strong inflows (and may even 
be outflows) observable at ~ kpc (Block et al. 2001; Garcia-Burillo 
et al. 2005; Haan et al. 2009; Stoklasova et al. 2009; Durbala et al. 
2009). Our simulations demonstrate that this can in fact be the case 
even in systems with large time-averaged accretion rates. 

The key parameter determining the qualitative behavior at 
each scale in our simulations is the ratio of the total mass in a ro- 
tationally supported disk (gas and stars) to the pre-existing mass in 
any spherical component (bulge, halo, or black hole); see Figure 13 
and § 5. There is, of course, a considerable literature studying the 
instability of self-gravitating disks in disk/bulge systems, which is 
applicable to the intermediate scales (~ 10— 100 pc) in our simu- 
lations (see e.g. Bardeen 1975; Athanassoula 2008; Narayan et al. 
1987; Raha et al. 1991; Christodoulou et al. 1995; Earn & Sell- 
wood 1995; Bournaud & Combes 2002; Dubinski et al. 2009). As 
outlined in Shlosman et al. (1989), most of these criteria amount 
to variations on the Ostriker-Peebles criterion that the kinetic en- 
ergy of rotation inside some radius be > 20% of the potential en- 
ergy. This requires increasing the surface density of disky material 



within the pre-existing bulge/halo effective radius to a value com- 
parable to, or equal to, that of the bulge/halo. For major, gas-rich 
mergers (/ gas > 0.3), this condition is almost always met. However, 
for gas-poor mergers, or for isolated barred galaxies, it is not clear 
how often systems will in fact trigger secondary instabilities. We 
find cases where weak large-scale bars trigger no significant subse- 
quent activity on smaller scales (Figure 6). 

In the gravitational potential of the central BH, the criterion 
for further inflow is somewhat altered, although we believe that it is 
qualitatively similar (§5.2): lopsided m = 1 modes are dynamically 
important provided that the total mass in the disky component is 
~ 0. 1 — 1 Mbh in the central ~ 10 parsecs (Figure 13; note the large 
scatter even at fixed disk to BH mass). We find that this criterion 
is generally satisfied, so long as the larger scales (~ 10— 100 pc) 
are themselves unstable. The critical barrier, then, to triggering a 
cascade of instabilities leading to significant BH accretion is the 
presence of a large self-gravitating gas inflow at ~ a kpc. 

There is significant evidence that galaxies in fact meet the cri- 
teria outlined here for AGN fueling. Studies of ongoing gas-rich 
mergers find that the gas densities reach and exceed the self-gravity 
threshold, with gas densities similar to those predicted here (Fig- 
ures 11-12; see e.g. Tacconi et al. 1999; Hibbard & Yun 1999; Laine 
et al. 2003; Iono et al. 2007; Schweizer & Seitzer 2007). Similarly, 
studies of early-type galaxies find that the stellar surface densities 
can be separated into a distinct "dissipational" component - the 
stars at small radii which must have formed in a dissipational star- 
burst - and a "dissipationless" envelope of violently relaxed stars 
(see Hopkins et al. 2008a). The starburst relic component is inferred 
to have formed from self-gravitating gas; indeed, this appears to 
define its size-mass relation (Hopkins et al. 2009a). In addition, 
the criterion for significant gas inflow given a modest pre-existing 
bulge corresponds to a surface density in the rotationally supported 
component of ~ 10 1 1 Mq kpc~ 2 at ~ lOpc (Figure 12), for an ~ L» 
system. This is comparable to the characteristic stellar surface den- 
sity observed at these radii in massive spheroids (see Hopkins et al. 
2009e, and references therein). 

One of the strong predictions of this work is that the link be- 
tween high gas surface densities and inflow leads to a correlation 
between BH accretion rate and star formation (§ 6). We have as- 
sumed a star formation law in which p, oc p 15 , which corresponds 
to a fixed star formation efficiency per local dynamical time. This 
is supported by models of turbulence-regulated star formation, in 
which the efficiency per dynamical time is a weak function of 
ambient conditions such as the Mach number of the turbulence 
(Krumholz & McKee 2005). With this assumption, we predict a 
correlation between BH accretion and the star formation at various 
scales in galactic nuclei, albeit one with significant scatter (Fig- 
ure 14). These predictions agree reasonably well with current ob- 
servations (Wang et al. 2007; Wu & Cao 2006; Silverman et al. 
2008). In the future, it will be interesting to explore the possibility 
that the star formation efficiency in galactic nuclei is significantly 
lower than the mean ISM value, as has been suggested by several 
authors (Thompson et al. 2005; Begelman & Shlosman 2009; Lar- 
son 2009). 

If the asymmetric stellar structures found here are long-lived 
(decay times of ~ Gyr), they should in principle be observable 
around nearby massive BHs. One tantalizing possibility is that the 
ubiquitous eccentric stellar and gaseous disks we find at ~ 1 — 10 
pc may explain the nuclear stellar disk seen in M31, which has 
been the subject of considerable study (see e.g. Lauer et al. 1993; 
Tremaine 1995; Salow & Statler 2001; Sambhus & Sridhar 2002; 
Peiris & Tremaine 2003; Bender et al. 2005). Hopkins & Quataert 
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(2010) examine this possibility in detail and show that the spa- 
tial scales ~ 1 — lOpc, moderate eccentricities, and stellar masses 
~ 0.1 — IMbh found in our simulations are all comparable to those 
observed. Eccentric nuclear disks have also been observed (albeit 
in less detail) around a number of other massive but inactive BHs, 
including NGC4486b (Lauer et al. 1996), M83 (Thatte et al. 2000), 
and VCC128 (Debattista et al. 2006). This strongly suggests that 
we can probe the the instabilities that produced the growth of mas- 
sive BHs in the "fossil" relics of that accretion. 

The dominant uncertainties in the models presented here are 
our treatment of the ISM physics, star formation, and feedback. 
Lacking a full micro-physical understanding of these processes - or 
the resolution to directly incorporate such a model - it is necessary 
to adopt a sub-resolution prescription for their impact. Our model 
is that stars form at all radii, in local, self-gravitating clumps, with 
a fixed efficiency relative to the local dynamical time. In addition, 
feedback from stars (supernovae, radiation, stellar winds/outflows, 
etc.) contributes, potentially along with cosmic rays and magnetic 
fields, to generating a turbulent velocity that is much larger than the 
thermal sound speed of the ISM. It is this turbulent "sound speed" 
that is included in our ISM equation of state (Figure 1). Observa- 
tions suggest that these assumptions are at least plausible even at 
the small scales we model here: a local Kennicutt-Schmidt relation 
holds in nuclear starbursts and on small scales in galaxies (Kenni- 
cutt 1998; Kennicutt et al. 2007; Bigiel et al. 2008) and star forma- 
tion efficiencies per dynamical time do not appear to evolve even in 
very dense regions (comparable to the highest gas densities mod- 
eled here; see, e.g., Tan et al. 2006; Krumholz & Tan 2007). More- 
over, AGN are observed to have associated nuclear star formation 
on scales as small as a pc; the densities of stars formed in situ at 
these scales (Davies et al. 2007; Hicks et al. 2009) are reasonably 
consistent with our predictions. Large turbulent and/or non-thermal 
gas velocity dispersions are also ubiquitous, from starburst nuclei 
on ~kpc scales to the sub-pc scales around AGN probed by water 
masers (Downes & Solomon 1998; Westmoquette et al. 2007; Tac- 
coni et al. 1999; Iono et al. 2007; Kondratko et al. 2005); these ob- 
servations motivate our choice of sub-resolution feedback models 
(Figure 1). To study the impact of ISM uncertainties on our results, 
we explicitly used different star formation efficiencies in some of 
our simulations; although this does have a quantitative affect on the 
resulting gas densities and inflow rates, the physical processes that 
we have identified as driving accretion remain the same. 

As noted in §2, the primary consequence of the large effective 
sound speed in our models is to increase the Jeans and Toomre 
masses, thus suppressing the formation of small-scale structure. 
That is, we effectively smooth over the smaller scales on which 
our treatment of the ISM physics is not appropriate. This implies 
that in our model, most of the gas, most of the time, resides in rel- 
atively diffuse structures, rather than being bound to very compact 
clusters, as would be the case if we did not include stellar feed- 
back and/or a sub-grid sound speed (see Appendix B). We suspect 
that this approach is reasonable for the global gravitational dynam- 
ics highlighted in this paper. In particular, our calculations indicate 
that over a remarkably large dynamic range, from 0. 1 pc to lOkpc, 
the torques and angular momentum transport are gravitational, and 
scale simply with the magnitude of the asymmetry in the potential, 
not with the sound speed of the gas (Fig. 9 and §4). In Paper II, we 
will present a more detailed comparison between analytic models 
and our numerical results that supports this conclusion. 

In our calculations, we have neglected AGN feedback in order 
to study how gas gets down to a BH in the first place. To the ex- 
tent that feedback is important on these scales, our simulations may 



better approximate the "early" stages of BH growth, when the BH 
is still relatively small (compared to, e.g., the still-forming bulge), 
and before it reaches a critical mass or luminosity at which feed- 
back becomes dynamically important. When BH growth becomes 
self-regulated, inflow and outflow are coupled, and the problem of 
AGN fueling must be addressed with a model that includes AGN 
feedback as a function of gas properties and spatial scale. 
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Figure Al. Example resolution test for a nuclear scale simulation. The 
same initial conditions re-simulated with the particle number shown (lower- 
resolution and nominal cases of Nf8h2blh; see Table 3). For > 0.5 X 10 
particles, our results are converged. 
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APPENDIX A: NUMERICAL TESTS 

We have performed a number of tests of our methodology, in an 
effort to ensure both that we have captured the most important 
physics in our simulations, and to ensure that there are no artifi- 
cial effects occurring due to numerical or resolution artifacts. We 
describe some of the key tests here. For general tests of GADGET- 
3, we refer the reader to Springel (2005); Springel & Hernquist 
(2003); Springel et al. (2005b); Springel & Hernquist (2002). 

Al Resolution Tests 

For the initial galaxy-scale simulations (Table 1), resolution tests 
have been extensively described previously (Di Matteo et al. 2005; 
Robertson et al. 2006b; Cox et al. 2006; Hopkins et al. 2009a); 
these range from ~ 10 5 particles to ~ 10 7 particles, and spatial res- 
olutions/softenings from ~ 20 pc to ~ lOOpc. For our purposes, 
the quantities of interest (e.g. the gas inflows to ~ lOOpc scales) 
are well-converged (see e.g. Appendix B of Hopkins et al. 2009a). 

We have performed a series of analogous resolution tests on 
our intermediate (Table 2) and small-scale simulations (Table 3), 
covering a similar span in particle number. For the intermediate- 
scale simulations, we have varied the SPH smoothing length and 
gravitational softening length from < 1 pc to ~ lOpc. We find that 
the results of particular importance for this paper, e.g., the total 
gas mass that loses angular momentum and flows to smaller scales, 
are well-converged for all the resolutions < lOpc used here. This 
is not surprising, because for physical reasons (discussed in the 
main text), the material reaches an angular momentum barrier at 
~ lOpc. For our small-scale simulations, we find similar behav- 
ior for all force resolutions <C 1 pc (experimenting with a range of 
softenings from 0.02 — 1 pc), and with good convergence at higher 
particle numbers > a few 10 5 . Given the range of resolution ex- 
plored, comparison with Figure 10 demonstrates that, especially 
in our high-resolution studies, even the verticle structure of the 



gaseous disks is well-resolved (in the highest-resolution cases, the 
typical h sm \/R < 0.01). Figure Al shows an example resolution test 
varying from 5 x 10 5 — 2 x 10 7 particles. 

A2 Do the Results Depend on Large-scale Tidal Fields? 

Because we are simulating the properties of large-scale non- 
axisymmetric systems, one might worry that the tidal forces from 
large radii could be important. This is certainly the case, e.g., in 
cosmological simulations. In the present context, the question is 
whether we need to include the matter distribution at ~ lOkpc in 
real-time to understand the dynamics of material at < lOOpc. 

We assess this directly as follows. For a given re-simulation 
of the nuclear regions, we determine the matter distribution of the 
larger-scale simulation from which the re-simulation is drawn and 
analytically fit the non-axisymmetric contributions to the potential 
at all radii. 10 Figure A2 compares the results from a re-simulation 
in which we include the large-scale non-axisymmetric potential at 
> kpc with our standard approach, in which we do not include this 
potential and in which the matter at large radii is not included in the 
re-simulation (so the information is lost). This is simulation If9b5 
in Table 2, a ~ 0.01 — 1 kpc re-simulation of the galactic nucleus 
during a galaxy-galaxy merger simulation at the peak of nuclear 
activity, just after the coalescence of the two galactic nuclei. At 
large scales, tidal tails are still present (these will take up to sev- 
eral Gyr to fully relax), and the system is largely unrelaxed outside 
of a few hundred pc, so this is a time when the large-scale effects 
would a priori be the most important. In the simulations shown in 
Figure A2 (left panel), we analytically fit the large-scale tidal field 
in the galaxy-galaxy merger simulation at multiple times, and in- 
terpolate it in time in our re-simulation. 

Figure A2 shows that there is effectively no difference be- 
tween the simulations with and without the large-scale tidal field, at 
any of the radii of interest. Including large-scale tidal forces does 
somewhat truncate the large spiral wave pattern that appears as a 
secondary instability in the intermediate-scale simulation, but this 
truncation occurs only on the largest radii > 500 pc, radii that we 
are not trying to model in detail in the re-simulation. Most impor- 
tantly, the scales which the re-simulation is intended to accurately 
represent, < lOOpc, are indistinguishable in the two simulations. 
The reason for this is ultimately that the mass profile is a relatively 
steep function of radius. If a mass SM is enclosed in an annulus at 
a radius r, then the tidal force it produces is oc |a|<5M/r 3 , where \a\ 
is the magnitude of the non-axisymmetric component of the mass 
distribution at that annulus. The contribution of a given large radius 
to the tidal force at small radii is thus oc p(f), the local mass den- 
sity. This, however, is always a decreasing function of (increasing) 
radius, so that material at larger radii makes little contribution to 
the tidal force. This can be compared with the role of tidal forces 
in determining the angular momentum of large-scale dark matter 
structures and ultimately (via these structures) dark matter halos; 
since the large-scale matter field has p ~constant on average, there 
are significant contributions to tidal forces from all radii. This is 
why including the information about large-scale tidal fields is crit- 
ical for simulations of cooling and galaxy formation, but it is not 
important when following the further inflows of material in indi- 
vidual galaxies. 



10 To do this, we adopt the radial basis expansion set proposed in Hernquist 
& Ostriker (1992), with standard spherical harmonics. This is similar to the 
density modeling described in § 2, and can be performed to arbitrary order, 
but only the first few terms are not noise-dominated. 
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Figure A2. Simulations testing the effects of large-scale tidal forces in our re-simulations. Left; Re-simulation of the central ~ lOpc-kpc of a larger-scale 
galaxy merger simulation, in which the full time-dependent potential from the large-scale simulation is analytically included in the re-simulation. We show the 
gas density and star formation rate density as in Figure 3, three large-scale dynamical times after the beginning of the simulation. Center: Same re-simulation, 
but without the additional potential from the large-scale material; i.e., the gravitational potential is only from the "live" material being re-simulated in the 
central kpc. The image is at the same point in time. Right: Gas density profiles of both re-simulations at the peak of the resulting inflows. Because the density 
declines rapidly with radius, the large-scale density/potential fields make no significant difference to the dynamics at the these scales. 



We have performed similar experiments to those shown here 
for other intermediate-scale simulations (both in mergers and iso- 
lated galaxies), and for our smallest-scale simulations (~ 0.1 — 
lOpc), and reach the same conclusions in each case. 

A3 Do the Results Depend on Initial Conditions? 

A reasonable question about our re-simulation technique is whether 
it might be sensitive to the precise initial conditions in the nuclear 
disk that is initialized at high resolution. In particular, there is lit- 
tle detailed information about the structure or dynamics of the nu- 
clear disk in the original larger-scale simulation, precisely because 
it lacked the necessary resolution to study the small-scale nuclear 
dynamics. Here we show that our results do not depend on the de- 
tails of how we initialize the re-simulations in the nuclear region. 
The fundamental reason for this is that the nuclear dynamics is gov- 
erned by instabilities that are self-consistently generated on small 
scales. These instabilities depend primarily on the gross structural 
properties of the nucleus, and rapidly lose memory of the initial 
conditions. 

Figure A3 shows the results of five different re-simulations, 
each with different initial conditions. These are all ~ lOOpc re- 
simulations of the central regions from a galaxy-merger simula- 
tion (intermediate-scale simulation If9b5 in Table 2), just after the 
coalescence of the two galaxy nuclei. Panel A is a rigorous re- 
simulation of the matter distribution in the central kpc of the large- 
scale simulation, with the matter initialized "as-is" from the large- 
scale simulation. Panel B is the same except we do not include 
any of the non-axisymmetric modes in the density, pressure, po- 
tential, etc. in the initial conditions; i.e., we azimuthally average 
the initial conditions used in Panel A. Panel C is even more ide- 
alized: the system is axisymmetric and with exponential disk and 
Hemquist profiles for the "disky" and "spherical" components in 
the re-simulation, respectively. The masses and effective radii of 
these components match the original simulation, but the profiles 
are these simple analytic approximations rather than the true pro- 
file from the larger-scale simulation. To test whether the presence of 
particular non-axisymmetries is important, Panel D includes non- 
axisymmetric structure in the initial conditions, but the precise 



modes are randomly chosen in phase and amplitude and so do not 
match the actual non-axisymmetries from Panel A. Finally, all of 
the above simulations include gas at small radii in the initial condi- 
tions; this is because the finite smoothing length in the larger-scale 
simulation spreads gas out over a region ~ 100 pc even though the 
dynamics below this scale is completely incorrect. In order to be as 
conservative as possible, Panel E considers an initial "hole" in the 
gas distribution at radii that are not resolved in the large scale simu- 
lation; for numerical reasons, the "hole" is a sharp power-law cutoff 
inside lOOpc. This could represent a very hard angular momentum 
barrier, such as a strong inner Linblad resonance. 

The images and surface density plots in Figure A3 are 10 8 yr 
after the beginning of the simulation. There is remarkably little dif- 
ference in the mode structure visible in the images. More quantita- 
tively, the surface density profiles are very similar, with significant 
gas inflow to ~ 10 pc. This demonstrates that the precise initial 
conditions used in the re-simulations are not important for most 
of our results. The fundamental reason for this is that the instabil- 
ities that dominate the dynamics arise from the internal structure 
of the material; they depend on the mass and scale-lengths of the 
rotationally and dispersion supported components, but not on the 
initial seed amplitudes of various modes. In the case of axisymmet- 
ric initial conditions (Panel B), the initial perturbations take slightly 
longer to appear because the initial power present to be amplified 
is smaller (it is due to particle noise); and in the case of strong ini- 
tial modes (Panel D), the instability is slightly more developed at 
the fixed time shown here. However, these time shifts have little 
long-term effect on the evolution of the system. 

The robustness of our conclusions to resolution, tidal fields, 
and initial conditions is also demonstrated by the fact that our ultra- 
high resolution simulations, discussed at length in the text, give 
identical results to our standard re-simulated runs. These ultra-high 
resolution simulations allow us to bridge the boundaries of our 
standard re-simulations, and so do not suffer the potential nega- 
tive effects considered here. Figure A4 shows an illustrative ex- 
ample of one such simulation about mid-way through the period 
of peak activity. The simulation is initialized as one of our typical 
intermediate-scale re-simulations, but with ~ 10 7 particles, giving 
a final resolution of < pc. On large scales, the dynamics is similar 



© 0000 RAS, MNRAS 000, 000-000 



How Do Massive Black Holes Get Their Gas ? 33 




Figure A3. Effects of initial conditions on our re-simulations of galactic nuclei. Panels A-E show the gas density and star formation rate density 10 8 yr after 
the beginning of the simulation. Bottom right: Surface density profiles of all five simulations at the same time. A: Re-simulation of the central ~kpc from a 
larger-scale merger simulation (as in Figure A2). The initial conditions are taken exactly from the larger-scale simulation, i.e., with inherited non-axisymmetric 
modes. B: Same as A but in this case we initialize azimuthally averaged (i.e., axisymmetric) profiles. C: Even more idealized: we initialize an axisymmetric 
bulge+disk+halo+BH system with analytic profile shapes; the mass and radius of each component is the same as that of the material in the original large-scale 
simulation at these radii, but the profile shapes are not. D: Same as A, but with random seed non-axisymmetric modes having order unity random amplitudes 
and phases E: Same as A, but the gas at < lOOpc is removed so that there is initially a central "hole" in the gas distribution. Comparison of A-E demonstrates 
that the precise initial conditions have little effect on the mode dynamics and gas transport. Morphologically, systems with stronger (weaker) initial seeds 
develop modes slightly more quickly (slowly), but this time shift has little significant or lasting effect on the dynamics. 



to our other intermediate-scale runs, with an m = 2 mode forming 
rapidly (there are non-trivial mode amplitudes from m = 1 — 3). 
The inflow leads to large gas masses at small radii; by the active 
phase shown in Figure A4, the lopsided, eccentric nuclear disk is 
evident. The properties of this disk are statistically indistinguish- 
able from those in our nuclear-scale re-simulations, given similar 
boundary conditions (compare e.g. Figure 2). 



APPENDIX B: ISM PHYSICS 

As discussed extensively in §2 & §7, the largest a priori uncer- 
tainty in our modeling is the question of what physics should be 
included to describe the behavior of the ISM. Our approach is to 
include a large "turbulent" sound speed in the equation state, as a 
subgrid model for the effects of feedback by stellar winds, radia- 
tion pressure, and supernovae on the ISM. We choose the turbulent 
velocity largely by comparison to observed systems (Fig. 1). The 
resulting turbulent velocities ~ 30— 100 kms -1 are also reason- 
ably consistent with theoretical estimates of the feedback required 
to maintain Q ~ 1 and disrupt molecular clouds in starburst envi- 
ronments (Thompson et al. 2005; Murray et al. 2009). 

To illustrate the importance of including a subgrid model, in 
this Appendix we compare the results of our standard models with 
a simulation in which the ISM is isothermal at ~ 10 4 K, i.e., c s ~ 
10 kms -1 (<5>eos = 0). We have experimented with a range of isother- 
mal cooling floors from c , = 1 — 15kms -1 (100 — 3 x 10 4 K), but 



the results here are generic to this range. Figure Bl shows sev- 
eral of the key results comparing the c, ~ 10 kms - ' simulation 
with our fiducial q cos = 0.125 simulation, for a re-simulation on 
~ lOOpc scales. From the images it is clear that after just a few lo- 
cal dynamical times in the intermediate-scale disk, the simulation 
without feedback has violently fragmented into large clumps. It is 
important to stress that in this simulation we fully resolve the Jeans 
and Toomre lengths so the fragmentation is physical, not numerical. 

The difficulty with the c s ~ 10 kms - (q ms — 0) simulation 
is that, absent a full model for how feedback can either suppress 
their formation or suppress their internal star formation, the only 
possibility is that the clumps "run away" and turn efficiently into 
stars. The runaway and final consequences of the fragmentation for 
star formation and the IMF are sensitive to details of the cooling 
rates and cooling function, but the end result of runaway star for- 
mation is inevitable if the gas is allowed to be arbitrarily cold and 
has a small cooling time (Gammie 2001; Nayakshin et al. 2007; 
Thompson et al. 2005). This has several dramatic effects. First, the 
global star formation efficiency is significantly enhanced - for the 
same mean gas surface density, most of the gas will turn into stars 
in ~ 1 — 4 local dynamical times absent feedback. This is a fac- 
tor of 10 — 50 higher than implied by the Schmidt-Kennicutt laws. 
For comparison, our fiducial q eos = 0. 125 simulation has a star for- 
mation efficiency per dynamical time of a few percent, much more 
consistent with observations. Because of the rapid and efficient star 
formation in the c s = 10 kms -1 simulation, the clumps quickly 
convert a large fraction of their mass to stars; the resulting mas- 
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Figure A4. Illustration of one of our ultra-high resolution simulations (Inf28b2h). The gas density distribution is shown face-on, as in Figure A3, in most 
panels, and edge-on at bottom right. The initial conditions are those of a typical intermediate-scale re-simulation, but with ~ 10 7 particles and sub-pc resolution. 
This obviates the need for a secondary "re-simulation" of the nuclear scales. The large-scale modes fuel gas inwards and here lead to a large nuclear gas mass, 
which then forms a lopsided m = 1 mode around the BH and drives further inflow. The eccentric nuclear gas disk is clearly visible, and similar to those in our 
nuclear-scale re-simulations; in fact there is no statistical difference, given similar boundary conditions. 



sive stellar clusters sink and coalesce via clump-clump scattering 
and/or dynamical friction, to the center of the galaxy. As the lower- 
right panel of Figure Bl shows, most of the mass flowing into the 
central ~ lOpc in the c s ~ 10 kms~' simulation is in the form of 
stars, rather than gas. This results in the formation of an extremely 



massive nuclear star cluster, with a mass of ~ 10 Mq inside the 
central ~ 10 — 20 pc - a factor of at least ~ 10 larger than ob- 
served nuclear star clusters (Boker et al. 2004; Cote et al. 2006; 
Ferrarese et al. 2006a). Indeed, the resulting stellar surface den- 
sity at small scales (left panel of Fig. Bl) significantly exceeds the 
highest nuclear densities observed in any cusp ellipticals (Ferrarese 
et al. 2006b; Kormendy et al. 2009; Lauer et al. 2007), or for that 
matter any nuclear star clusters, globular clusters, nuclear disks, or 
other high-density stellar systems (Hopkins et al. 2009e). 

In contrast to the constant c s ~ 10 kms - simulation, our fidu- 
cial simulation with significant subgrid feedback does not clump 
up into many dense star clusters. This is also not completely phys- 
ical since in reality we expect that the ISM should have signifi- 
cant sub-structure on the scales we model (including, e.g., molec- 
ular clouds and star clusters). However, when such self-gravitating 
gaseous clumps form, observations strongly suggest that are likely 
to be short-lived and/or to inefficiently form stars. Our subgrid 
model assumes that most of the mass remains in a diffuse ISM, 
rather than becoming bound in dense star clusters. Not only is the 



latter physically implausible, but Figure Bl demonstrates that it 
strongly violates observational constraints. Note that there are con- 
ditions under which even an extremely cold gaseous disk could 
avoid catastrophic fragmentation (see e.g. Lodato & Rice 2004; 
Rice et al. 2005; Boley et al. 2006; Krumholz et al. 2007), partic- 
ularly the case in which the cooling time is sufficiently long such 
that continued cooling is offset by gravitational heating (giving rise 
to sustained local structure such as tightly-wound spiral arms, but 
not runaway fragmentation). However, the cooling rates under typ- 
ical conditions in our simulations are much too rapid to reach this 
regime; we could, of course, modify the cooling function to sup- 
press the cooling rate and/or mimic some continuous heating term 
- but this is effectively equivalent to adopting a new sub-grid feed- 
back/microphysics model, and we find it has the same qualitative 
effects. 

Figure B2 shows the effects of varying the efficiency of feed- 
back from supernovae and massive stars, encapsulated in the pa- 
rameter ^cos- Specifically, we show results of a set of nuclear scale 
simulations (Nf8h2b4q in Table 3) having q tos = 0.02 - 0.40 (ef- 
fective turbulent velocities ~ 20— 100km s _1 ), roughly the lower 
and upper limits allowed by observational constraints for the sys- 
tems of interest in Figure 1. Figure B2 shows that the choice of 
sub-resolution model has a significant effect on the amount of re- 
solved sub-structure in the simulation. This is not surprising since 
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Figure Bl. Examples of the problems that occur if systems are evolved without some "feedback" to prevent runaway cooling and fragmentation. Top Left: 
Gas density in such an initially smooth re-simulation of the central ~kpc (unlike in the main text, here darker areas represent higher-density to highlight the 
fragmentation). The system here can cool to 10 4 K and has no feedback (q eos = 0); the image is shown after ~2x 10 6 years (less than the global dynamical 
time). Top Right: The same initial conditions but for our fiducial model - a moderate subgrid ("turbulent") sound speed ~ 30kms — 1 (q e os = 0.125) in dense 
star-forming regions. In the no feedback g cos = case, gas clumps at the leans scale and rapidly turns into stars - leading to a severe violation of the observed 
global Kennicutt (1998) relation. Bottom Left: Stellar density profiles after 10 7 yr, for simulations with no feedback (two examples shown, with different 
cooling floors as labeled) and our fiducial g cos = 0.125 model; the initial density is just 10 10 Mq kpc -2 at small radii. We also show the observed stellar mass 
density profiles of massive cusp ellipticals in Virgo from Kormendy et al. (2009), chosen to have the same mass at > 0.5 — lkpc. In the absence of feedback, 
runaway fragmentation and the inability of clumps to dissolve inevitably leads to the formation of an extremely massive nuclear stellar cluster with a mass 
and surface density at least an order of magnitude larger than any observed. Bottom Right: The accretion rate into the central lOpc for both gas (dotted) and 
already-formed stars (solid). In the q eos = models, sinking clumps provide large gas accretion rates ~ 10 — 100 Mq yr — 1 , but the stellar inflow rate typically 
exceeds the gas inflow rate by a factor of ~ 10; this is inconsistent with the observed stellar densities at small radii. Our fiducial, moderate-feedback case, by 
contrast, drives primarily gas - not stars - to small radii. 



larger turbulent velocities raise the Jeans mass/length, above which 
gravity is the dominant force. The key point, however, is that all of 
the simulations show a similar nuclear lopsided disk. In terms of the 
properties of inflow and global instability discussed throughout this 
paper, the differences produced by changing q tos are similar to the 
differences produced by somewhat different galaxy properties. The 
fundamental reason for the weak dependence on the subgrid model 
is that the torques in our simulations are primarily determined by 
gravity, not hydrodynamic forces or viscosity. The primary role of 
the subgrid feedback model is simply to prevent catastrophic frag- 
mentation of the galactic gas. Provided this can be accomplished, 
our qualitative conclusions do not depend on the details of the feed- 
back model. 
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Figure B2. Effects of varying the assumed efficiency of stellar feedback on the resulting nuclear eccentric disk, given identical inflow conditions at lOOpc. 
We show the gas surface density with color indicating specific SFR (as Figure A3), for several choices of q ms (from the survey Nf8h2b4q). The range shown 
corresponds to sub-grid turbulent velocities ranging from ~ 20 — 100 km s in the high star-formation rate regions of the ISM. The efficiency of feedback 
affects the level of substructure. However, the general character and appearance of the lopsided disk mode is similar over the entire plausible range. Provided 
the catastrophic fragmentation in Figure B 1 can be avoided, our qualitative conclusions are independent of g cos . 
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